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Abstract 

Phase field models for two-phase flow with a surfactant soluble in 
possibly both fluids are derived from balance equations and an energy 
inequality so that thermodynamic consistency is guaranteed. Via a formal 
asymptotic analysis, they are related to sharp interface models. Both cases 
of dynamic as well as instantaneous adsorption are covered. Flexibility 
with respect to the choice of bulk and surface free energies allows to 
realise various isotherms and relations of state between surface tension 
and surfactant. Some numerical simulations display the effectiveness of 
the presented approach. 

1 Introduction 

Surface active agents (surfactants) reduce the surface tension of fluid interfaces 
and, via surface tension gradients, can lead to tangential forces resulting in the 
Marangoni effect. Biological systems take advantage of their impact on fluids 
with interfaces, but surfactants are also important for industrial applications 
such as processes of emulsification or mixing. While often much experience 
and knowledge is available on how surfactants influence the rheology of multi- 
phase fluids, the goal is to understand how exactly the presence of a surfactant 
influences coalescence and segregation of droplets. 

Surfactants can be soluble in at least one of the fluid phases and the exchange 
of surfactants between the bulk phases and the fluid interfaces is governed by 
the process of adsorption and desorption. Ward and Tordai (54] derived a time- 
dependent relation for the surfactant density at the interface and the surfactant 
density at the adjacent bulk phase (known as the sub- layer or sub-surface). To 
compute the interfacial density, a closure relation between the two quantities has 
been proposed in the form of several different equilibrium isotherms [18] 1331 132j , 
where the underlying assumption is that the interface is in equilibrium with the 
sub-layer at all times. This corresponds to the case of diffusion-limited adsorp- 
tion studied in Diamant and Andelman [16] . where the process of adsorption 
to the interface is fast compared to the kinetics in the bulk phases. However, 
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instantaneous adsorption is not valid in the context of ionic surfactant systems 
[IB] or when the diffusion is not limited to a thin layer [T^J EH H] • Therefore, 
we would like to be able to account for non-instantaneous adsorption in our 
models. 

Two-phase flow with surfactant is classically modelled with moving hyper- 
surfaces describing the interfaces separating the two fluids. We will derive the 
following sharp interface model for a domain O containing two fluids of different 
mass densities. We denote by ^^'(t), ^ 2 '(t) the domains of the fluids which 
are separated by an interface T(t): 

V-t> = inO (l) (i), (f.f) 
dt{p {i) v) + W-(p {l) v<E>v)^W-^- P I + 2ri^D{v)^ in W (i), (1.2) 

9*c«=V-(M«VG^(c«)) infiW(i), (1.3) 

[v]f = 0, v-u = u T on T(t), (1.4) 

\ P I-2ri {l) D(v)]lv = cj{c T )KV + V V (j(c T ) on T(i), (1.5) 

9 t , c r + c r Vr-w-Vr-(il/ r Vr7'(c r )) = [A/«VG' t (c«)]2iy on T(i), (1.6) 

a«(-l) i M«V^(c«)-iv = -( 7 '(c r )-^(c«)) onr(t). (1.7) 

Here v denotes the fluid velocity, is the constant mass density for fluid i, 
17 W is the viscosity of fluid i, D(v) = ^(Vu + ^r) 1 ) is the rate of deformation 
tensor, p is the pressure, I is the identity tensor, (?*(•) = d t (-) + v ■ V(-) is the 
material derivative, is the bulk density of surfactant in fluid i, Mj z ' is the 
mobility of surfactants in fluid i, Gj(cW) is the bulk free energy density associ- 
ated to the bulk surfactant in fluid i. On the interface, ur is the normal velocity, 
v is the unit normal on T pointing into Q( 2 \ c r is the interfacial surfactant den- 
sity, cr(c r ) is the density dependent surface tension, k is the mean curvature 
of r, Vr is the surface gradient operator, Vr - is the surface divergence, Mr 
is the mobility of the interfacial surfactants, 7(c r ) is the free energy density 
associated to the interfacial surfactant, and cr^ > is a kinetic factor that re- 
lates to the speed of adsorption. The above model satisfies the second law of 
thermodynamics in an isothermal situation in the form of an energy dissipation 
inequality. 

Equations (jl.ip and (|1.2p are the classical incomprcssibility condition and 
momentum equation, respectively. The mass balance equation for bulk surfac- 
tants is given by (jl.3p . Equation (|1.4j) states that the interface is transported 
with the flow and that not only the normal components but also the tangential 
components of the velocity field match up. The force balance on the interface 
(|1.5|) relates the jump in the stress tensor across the interface to the surface 
tension force and the Marangoni force at the interface. The mass balance of the 
interfacial surfactants is given by (|1.6[) . and the closure condition (|1.7[) tells us 
whether adsorption is instantaneous (aW=0, an isotherm is obtained) or dy- 
namic (cjW > 0, the mass flux into the interface is proportional to the difference 
in chemical potentials). 

The model studied in [TU] bears the most resemblance to the above model, 
where the setting of these papers is the diffusion-limited regime with a surfactant 
which is soluble in one phase only and (|1.7p is replaced by the relation 



7 '(c r HG'(c) c r = .g(c):=(7')- 1 (G'(c)), (1.8) 
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in which g plays the role of the equilibrium isotherm and where G is the bulk free 
energy of the phase in which the surfactant is soluble. Our approach is based 
on a free energy formulation, originated from [161 117] . where we gain access 
to equilibrium isotherms by setting = and choosing suitable functions for 
7 and Gi. Furthermore, for positive values of a'*' we are able to include the 
dynamics of non-equilibrium adsorption. 

The governing equations (jl.ip — (|1.7p form a free boundary problem. The 
phase boundary T(t) is unknown a priori and hence must be computed as part of 
the solution. Much previous work have been dedicated to explicitly tracking and 
capturing the interface using various numerical methods [57j [28] [56] [36] [43] [30] . 
However, the sharp interface description breaks down when topological changes 
occur. Phenomena such as breakup of fluid droplets, reconnection of fluid inter- 
faces and tip-streaming driven by Marangoni forces [2TJ [35J EH] involve changes 
in the topology of the interface. Numerically, complications also arise when the 
shape of the interface becomes complicated or exhibits self- intersections. These 
difficulties have led to the development of diffuse interface or phase field models 
to provide an alternative description of fluid/fluid interfaces. 

At the core of these models, the sharp interface is replaced by an interfacial 
layer of finite width and an order parameter is used to distinguish between the 
bulk fluids and interfacial layer. The order parameter takes distinct constant 
values in each of the bulk fluids and varies smoothly across the narrow interfacial 
layer. The original sharp interface can then be represented as the zero level set 
of the order parameter, thus allowing different level sets to exhibit different 
topologies. 

The width of the interfacial layer is characterised by the length scale over 
which the order parameter varies from its values at the bulk regions. The phase 
field model can be related to the sharp interface model in the asymptotic limit 
in which this width is small compared to the length scales associated to the bulk 
regions. Hence one can also view the phase field methodology purely as a tool for 
approximating the sharp interface equations. If the objective is to ensure that, 
in the limit of vanishing interfacial thickness, certain sharp interface models are 
recovered then there is a lot of freedom in constructing phase field models to 
meet one's needs (see e.g. EUJ). 

The review [1] provides an overview on diffuse interface methods in the 
context of fluid flows. In [55] H7J it was already proposed to combine a Cahn- 
Hilliard equation for distinguishing the two phases with a Navier-Stokes system. 
An additional term was included in the momentum equation to model the sur- 
face contributions to forces. In the case of different densities, Lowengrub and 
Truskinovsky [5T] derived quasi-incompressible models, where the fluid velocity 
is not divergence free. On the other hand, Abels, Garcke and Grim [T] derived 
a thermodynamically consistent diffuse interface model for two-phase flow with 
different densities and with solenoidal fluid velocities. Following the derivation 
in pp, we will derive three diffuse interface models, which approximate the sharp 
interface models in the diffuse-limited regime. 

More precisely, for the case of non-instantaneous adsorption (a^' >0), we 
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will derive the following model (denoted Model A) 

V-« = 0, (1.9) 
d t {pv)+V- {pv®v)=V- (-pI + 2r]{ip)D(v) + v<g) P m ~P (1) m(y)V/ij (1.10) 
+ V ■{Ka{c T ){5{p,Vp)I-eVp®Vp)), 

a;<p=v-(m(<p)Vn), (l.n) 

p + V-(Kea(c r )Vp) = -<j(c r )W\p) + £ £^)(Gi(c«) -GKc W )c (<) ), 
e * — ' 

i=1.2 

(1.12) 

a; & (^) c w ) = v • (m c w ( c « vg$ (c w )) (1.13) 

+ /3W^,V^)( 7 '( C r )-GUc«)), 1 = 1,2, 
d'{K5{p,Vp)c T )^V ■ (j\I v {c T )K5{p,V p)Vi {c T )} (1.14) 

-%,V^)^/3W( 7 '(c r )-G:( C W)). 

j=l,2 

Here £ is a length scale associated with the interfacial width, cp is the order 
parameter that distinguishes the two bulk phases. In fact ip takes values close 
to ±1 in the two phases and rapidly changes from —1 to 1 in an interfacial 
layer. The functions £i(ip) and 8{p,Vp) act as regularisation to the indicator 
functions of fiW and T, respectively. The quantity ^ l '=K/w- % ' is related to 
the adsorption kinetics and A is a constant. Equations (|1.9[) and (jl.lpp are the 
incomprcssibility condition and the phase field momentum equations, respec- 
tively. Equation (jl.llj) together with (|1.12l) governs how the order parameter 
evolves and equations (|1.13[) and (|1.14p are the bulk and interfacial surfactant 
equations, respectively. 

We derive two additional models for instantaneous adsorption (a" =0): 
Model B models the case where the surfactant is soluble in only one of the 
bulk phases. It consists of (|1.9p — (|1.12p and replaces the bulk and interface 
surfactant equations (|1.13p , Q1.14p with 

3t{£{ip)c+K5g(c))- V-(M(c)^)VG'(c))- V- {M T {g{c))K8VG'{c)) = Q, 

(1.15) 

where g(c) is the adsorption relation between interface and bulk densities as in 

(USD- " 

The case where the surfactant is soluble in both bulk phases is covered by 
Model C, which consists of (fTUJ) - (fTTT2"|) and 

d' t (f i ( V)c (1) (q) + 6 {p)c {2) (q) + KSc r (q)) - £ V • (Mj (c« (q) % (p) Vq) 

i=l,2 

(1.16) 

-V-(M r {c r (q))KSVq) = 0. 

Here, q denotes a chemical potential where, as will be discussed in Section 3, 
we can express the surfactant densities as functions of q. 

The Model A is related to the approach in [35] . We modify the approach of 
|48| in such a way that an energy inequality is valid and such that we recover 
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the isotherm relations for adsorption phenomena in the limit of instantaneous 
adsorption. We deepen the asymptotic analysis in that it works with the orig- 
inal equation for the surface quantity and does not require the assumption of 
extending the surface quantity continuously in normal direction. Phase field 
models of surfactant adsorption that utilise the free energy approach of [16j [17] 
can be traced back to the models of jHH \EH H3] , where the latter is extended 
in [3U] and solved using lattice Boltzmann methods. The issue of ill-posedncss 
of the model is discussed in [3U] and three alternatives have been suggested. 
Phase field models that look into the behaviour of equilibrium configurations of 
fluid-surfactant systems can be found in (23] [50J and a detailed comparison of 
previous phase field models can be found in [38]. 

The structure of this article is as follows: In Section 2 we will derive the 
sharp interface model — (|1.7[) from basic conservation laws. We show that 
the sharp interface model satisfies a local energy inequality and present the func- 
tional forms for 7 and G that lead to five of the popular adsorption isotherms 
when = 0, namely those of Henry, Langmuir, Volmer, Frumkin and Fre- 
undlich. In Section 3, we present the derivation of phase field models based on 
the Lagrange multiplier method presented in [T] and show all of them satisfy a 
local dissipation inequality In Section 4 we show, via formally matched asymp- 
totics, that we recover (jl.ip — p.7p from Model A and (|1.8p from Models B and 
C in the limit e~ s-0. In addition, Model A can be shown to converge to the 
sharp interface problem with instantaneous adsorption when the kinetic term is 
chosen appropriately. In Section 5, we present ID and 2D numerics to support 
the asymptotic analysis. 

2 Sharp interface model 
2.1 Balance equations 

We consider a domain fl C K d , d = 1 , 2, 3, containing two immiscible, incompress- 
ible Newtonian fluids with possibly different constant mass densities p ,£ = 1,2. 
The domain occupied by the fluid with density is labelled as fl^ cKxIR d 
where we set (t) := {x e il; (t,x) G il 1 -^}. The two domains are separated by 
an interface T which is a hypersurface in M x R d such that r(i)ndf2 = where 
r(<) := {x £ O; (t,x) € T}. A surfactant is present which alters the surface tension 
by adsorbing to the fluid interface and, provided it is soluble in the correspond- 
ing fluid, it is subject to diffusion in the phases fl^\ We denote the fluid velocity 
field by v, the pressure by p, the bulk surfactant densities by c^,i= 1,2, and 
the interface surfactant density by c r . 

Balance of mass and linear momentum inside the phases lead to the following 
equations 

V ■ v = 0, ft (7j (i) v) + V ■ (7i w v <g> v) = d' t (p {l) v) = V • T (l) , 

where d' denotes the material derivative and T^>, £=1,2, is the symmetric 
stress tensor (due to conservation of angular momentum) . These equation hold 
in Qy-'{t)UQS 2 '{t). We assume that the two fluids do not undergo phase transi- 
tions and the phase boundary T(t) is purely transported with the flow where we 
also assume that there is no-slip at the interface, hence the tangential velocities 
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match: 

[u]i=0, v-i/ = ur- 

Here [-ff denotes the jump of the quantity in brackets across T from f^ 1 ) to 
v is the unit outward normal of T(t) pointing into Q( 2 \t), and Mr is the 
normal velocity of the interface. 

Let V(t) be an arbitrary material test volume in f2 with external unit normal 
u ext of V(t)nf2. If V(t)P\T(t) is non-empty then we denote its external unit 
co-normal by fi and write v^ t for the external unit normal of V (t) flfi^ (t) , 
i = l,2. In the bulk fluid regions, surfactants will be subjected to transport 
mechanisms consisting of only diffusion and convection. Hence, mass balance 
for bulk surfactants in a material test volume V(t) away from the interface T(t) 
yields 



d_ 

df 



c« = -/ J^-v ext 

v(t) Jav(t) 



where Jc is the molecular flux. By Reynold's transport theorem and using 
that V • v = 0, this leads to the pointwise law 

a t *cW+V-jW=0, i = l,2. (2.1) 

For a test volume V(t) intersecting T(t), we postulate 

c®+[ c r | 

v(t)nnw(t) Jv(t)nv(t) 




8(v(t)nn«(t))\r(t) ./a(v(t)nr(t)) 

where «7r is the interfacial molecular flux, tangential to T. Using Reynold's 
transport theorem, the surface transport theorem and the surface divergence 
theorem (see [B]) we obtain 




v(t)nnw(t) Jv{t)nv(t) 



d'c^+f (d'c r + c r V r -v) 

Jv(t)nr(t) 



'v(*)n!)W(t) J^(f)nr(t 
for the left hand side and 

i=i,2 •'8(v(t)nnw)(i))\r(t) ./8(v(i)nr(t)) 

1=12 Jd(v(t)rfoW(t)) Jv(t)nr(t) 

for the right hand side. Hence, using (|2.1[) the mass balance (|2.2p yields the 
following pointwise law for the interfacial surfactant: 

d'c r + c r V r • v = - V r • Jr + <L4D, 
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where qad = — [Jc^]\v is the mass flux for the transfer of surfactant to the 
interface from the adjacent sub-layers. When the mass flux qad is zero and the 
interfacial molecular flux is modelled by Fick's law, Jr = — Z? s Vrc r , we obtain 
the mass balance equation in |55j . 

2.2 Energy inequality 

We postulate a total energy of the form 

f [^\v\ 2 + G 1 (cW)}+[ [^ M 2 + G 2 (c< 2 >)]+/ 7 (c r ), (2.3) 
Jo(D(t) Jn( 2 )(t) Jr(t) 

where G\,G2 are the bulk free energy densities, and 7 is a surface free energy 
density. We assume that 7" > and G" > 0. The Legendre transform of the 
surface energy density then is well defined, and the density dependent surface 
tension u(c r ) is defined as 

a( C r ):=7(c r )-c r 7'(c r ). (2.4) 
Let V(t) be an arbitrary material test volume. Then 

|(E/ C^\v\ 2 + G^))+f 7 (c r )) 

at \ 7Ti Jv(t)nnW(i) Jv(t)nr(t) I 



=E/ (p®vd' t v + G' i (cU)d; C W) + f ( 7 '( c r )9 t 'c r + 7 ( C r )Vr-«) 

i=1 ./vnn« v ' Jvr\Y 

2 „ 

= J2 / fv.((TW) i D-G;(c (,) )Jc t ' ) )-2 l(,) : Vv + VG' 4 (c W ) • J«) 

J=1 iv(t)nn(')(t) v ' 

+ [ ~f'(c r )(-V r -Jr+qAD)+<j(c r )V r -v 
Jv(t)nr(t) 

2 

= V [ -T {i) : + VG ■ (c W ) ■ J® + [ - 7 '(c r )Jr-/x 

i=1 Jv(t)nci(')(t) Jd(y(t)nr(t)) 

+ £ / ((TW^v-G'M^jU).^ 

+ / ((TW)^-G' 1 (c( 1 »)jW). I , + ((T( 2 ))^-Gi(c( 2 ))j( 2 )).(-^ 

Jv(t)nr{t) 

+ [ Jr-Vr7'(c r )+ 7 '(c r )(Ji 1) -iv-j( 2 )-iy)+a(c r )Vr--y. 

Jv(t)nr(t) 

Decomposing the velocity field v on T(t) into its normal and tangential compo- 
nents 

v = uyv + v t , 

then gives 

/ tr(c r )Vr • (uyv + v t ) = / cr(c r )(Vrwr • I' + urVr • t'+Vr -w T ) 
iy(t)nr(t) iv(t)nr(t) s » ' y v ' 

—0 — nur 

= / — a(c r )nur — Vr&(c r )-v + / a(c r )f T -/2, 

./vmnrm ia(y(t)nr(t)) 
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where k = — Vr • v is the mean curvature and we have used integration by parts 
to obtain the last equality. Altogether we have 



£(£/ [^H 2 + G 4 ( C «)]+/ 7 (c r ) 

at \i=i Jv(t)na^(t) Jv(t)nT(t) ) 

E/ ((t«)^-g;(cW)j«)-^ 

^•/<9(v(t)nnw(t))\r(t) 

+ / (-7'(c r )J r -/x + cr(c r )i; T -/x) 

J9(v(t)nr(t)) 

2 . 

+ £ / f-rW: V-u + VGKc W )-Ji l) ) 

+ / J r -V r7 '(c r ) 

JV(t)nr(t) 

Jv(t)nr(t) v 7 
+ / (T (1) /y-«-T (2) iy-«-cr(c r )Kt;-iy-Vrcr(c r )-v) . 



>v(t)nr(t) 
Hence, if 

jW-VG^(c (i) )<0, innW(t), i = l,2, 

T (i) :Vv>0, inO (i) (t), i = 1,2, 

J r -Vr7'(c r )<0, on r(t), 

(J< 1 )- I ,)( 7 '(c r )-G' 1 ( C (1) ))<0, onr(t), 

(-j( 2 )^)( 7 '(c r )-G' 2 ( C ( 2 )))<0, onr(t), 

(-[T] 2 ^-cr(c r )Kiv-Vrcr(c r ))-u<0, on 

then we obtain the following energy inequality: 

!(£/ (^M 2 + G l(c «)) + / 7 (c r )) 

at \ l=1 iy(t)no< i )(t) Jv(t)nr(t) / 



l=1 \Jd(v(t)nn(*)(t))\r(t) , 

+ (- 1 '(c r )J r -fi + ( j(c r )v T -fi), 

Jd(v(t)nr(t)) 



id(v(t)nr(t)) 

where the right hand side represents the working on the arbitrary material test 
volume V(t) and the inequality indicates that the dissipation is non-positive, 
thus guaranteeing thermodynamic consistency |24l I26j . 
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2.3 General model 

Wc make the following constitutive assumptions: 

jW = -AfW(cW)VG£(c<*>), 
Jr = -Mr(c r )Vr7'(c r ), 

a^{c v ,S){-l) i+l J^-u = -{ 1 \c v )-G' i {S)), (2.5) 
T (i) = -pI + 2r 1 (l) D{v) 1 
-[T}lh> = <T{c v )ni> + Vr<j{c r ), 

where M c (4) (c«) > 0, A/ r (c r )>0, and a^{cF,c^)>0. 

The formulation presented in (|2.5|) utilises a free energy approach, first ap- 
plied to the kinetics of surfactant adsorption in |16l I17| , to model instantaneous 
adsorption kinetics. At adsorption/desorption equilibrium, the chemical poten- 
tials 7'(c r ) and G'(c) must be equal [551 001 [53] an d thus this approach allows 
us to cover the adsorption isotherms often used in the literature by selecting 
suitable functional forms for 7 and G. Hence, a"' > can be seen as a kinetic 
factor which relates the speed of adsorption to the interface or desorption from 
the interface to the deviation from local thermodynamical equilibrium. Let us 
summarise the governing equations of the general model for two-phase flow with 
soluble surfactant: 

Balance equations in fiW(i), i = l,2 : 

V-u = 0, (2.6) 
d t {p {i) v) + V • (pi- 2ij®D(v) +p {l) v<E> v) = 0, (2.7) 
d t *c w - V • (AfWVGj(cW)) = 0. (2.8) 

Free boundary conditions on T(t): 

[v]i = 0, v-v = ur, (2.9) 

[pjfu -2[r) (i) D(v)]\v = o{c t )kv + V r cr(c r ), (2.10) 

d' t c T + c r V r -v = Vr • (M r Vr 7 '(c r )) + [M« VG^(c«)]?^, (2.11) 

a (j) (-l) l Mi l) VG^(c W )-i/ = -(7 , (c r )-G^(c W )). (2.12) 

In this model, the surface tension cr : M + — ► M + is a (usually decreasing) func- 
tion of the surfactant density c r . The phenomenon known as Marangoni effect, 
where tangential stress at the phase boundary leads to flows along the interface, 
is incorporated into the model via the surface gradient of a in the momentum 
jump free boundary condition. 

2.4 Specific models 
2.4.1 Fick's law for fluxes 

By appropriate choice of the mobilities we obtain Fick's law for the surfactant 
both in the bulk and on the surface. If we set 

M«(cW) = i)W — * M T {c v ) = D T - 



G^(cW)' iV ; 7"(c r )' 
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for constant Fickian diffusivities Dc ,-Dr > 0. Then 

J® = -D®Vc®, J r = -D r V r c r . 

2.4.2 Instantaneous adsorption and local equilibrium 

We may assume that the process of adsorption of surfactant at the interface is 
instantaneous, i.e. fast compared to the timescale of convectivc and diffusive 
transport. This local equilibrium corresponds to the case that the bulk chemical 
potential G"(c) and the interface chemical potential 7'(c r ) are equal, i.e. we set 
a = in (|2.5p (we here only consider one of the bulk phases adjacent to the 
interface and, for simplicity, drop the upper index (i)). We obtain the following 
relation (also see [51 HOjl : 

7 '(c r HG'( C ) <=s> c T = g{c):={ 1 ')-\G'{c)), (2.13) 

where g : R+ — > R+ is strictly increasing. This function g plays the role of various 
adsorption isotherms which state the equilibrium relations between the two 
densities. 

Table [T] displays the functional forms for 7 and G in order to obtain the 
adsorption isotherms of Henry, Langmuir, Freundlich, Volmcr and Frumkin. 
The free energies are (variants of) ideal solutions. Here, c£ f is the maximum 
surfactant density on the interface, K a constant relating the surface density 
to the bulk density in equilibrium, <jq denotes the surface tension of a clean 
interface, B essentially is the sensitivity of the surface tension to surfactant, A 
in the Frumkin isotherm is known as surface interaction parameter while, in the 
Freundlich isotherm, A c measures the adsorbent capacity and N is the intensity 
of adsorption. 

2.4.3 Insoluble surfactants 

Neglecting (|2.8|) . ({2.120 and the jump term in ({2. lip gives a two-phase flow 
model with insoluble surfactant. 

2.5 Reformulation of the surfactant equations 

The strong form of the surfactant equations (|2.8|) . (|2.11|) . (|2.12|l can be reformu- 
lated into an equivalent distributional form using a result from Alt [3] . Let Xnw 
and Sr denote the distributions given by the indicator functions on and T 
respectively, see the Appendix for a precise definition. We now define 

h = ^(-y'(c r )-G[(cM)), J2 = -L( 7 '( c r)- G ^( c ( 2 ))). 

In the Appendix we show that 

dt (xnw ) + V • Otow C W« - xnw M« VG' X (c« )) = S rjl , (2.14) 
dt (xn< 2) ) + V • (xnw c^v- Xam M™VG' 2 (c^j) = S rj2 , (2.15) 
d t (S r c r ) + V • (6 r c r v - M r 5 r V 7 '(c r )) = -8 r (ji +32), (2.16) 
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a-cr Q 


-Bc T 


Bc^log(l-^) 


Isotherm 


Freundlich 


Volmer 




Relation 


/ r \ N 


ifc = r cr exp ( r c r ) 


7(c r )-(T 


«Bc r (logi-l) 


BcT log 


G(c) 


Bc(hg(A* Kc) - 1) 


Bclog(Xc) 


a-cr Q 




r r 

C M C 



Isotherm 


Frumkin 


Relation 


^c- c ,;_ cr exp( % ) 


7(c r )-(T 




flog^+c^log(l-f ))-^ 


G(c) 


Bc(log(Kc)-l) 


a-a 


^ + ^ M log(l-g 



Table 1: Possible functional forms for 7 and G to obtain the most frcque 
used adsorption isotherms and equations of state. 
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interpreted in its distributional formulation are equivalent to 

9 t c (1) +V ■(c^v-M^WG' 1 (c^)) = 0, inftW, 

AfWVGi(cW)-v = ji, on T, 
d t c^+V-(c^v-M^VG' 2 (c^)) = 0, inf2( 2 ), 
-M^ 2) VG' 2 {c^)-v^j 2 , onT 

and (|2.1ip respectively. 



2.6 Non-dimensional evolution equations 

To derive equations in a dimensionless form we pick a length scale L, a time scale 
T (or, equivalently, a scale for the velocity V = L/T), a scale E for the surface 
tension, and let C r = L~ 2 ,C = L~ 3 denote scales for the surfactant densities in 
the interface and in the bulk, respectively. 

The Reynolds number, as the ratio of advective to viscous forces, is defined 
as Re:= (~p^L 2 )/(rj^T). The capillary number, as the ratio of viscous to sur- 
face tension forces, is defined as Ca— (r]( 2 >L)/(TY,). Scaling the pressure by 
T 2 / L 2 ) we arrive at the following dimensionless fluid equations 



V*-t;* = 0, 



2f? d 



d t , (p v ») + V* • p*I— — — D(v*) + p v*®v 



Re 



p.I-^-D(v«) 
Re 



= 0, 



[i>*] 2 = 0, v*-v = u T „, 



1 



ReCa 



(tr*Ki/+Vr»o'*), 



where 17+ = 1, rf~ = rf^ /rf 2 ^ , p + = l, 75~ =p (1) /p (2) - Let 
7 



7* 



G« * = -4- , M c (i i = M„ w STL 3 , M r * = M r STL 2 



(2.17) 
(2.18) 
(2.19) 
(2.20) 



where 7*,Gi i * denote the dimensionless free energies and M^Mr,* denote the 
dimensionless mobilities. The dimensionless surfactant equations are given by 



9 t * cS° - V* • (M#V.G< ) =0 



r (i) 



CO- 



9* <£+<£V r . V r . ■ (Mr,*V r .7i(<£)) = M$V*G^( C W) 



(2.21) 
(2.22) 



aW(-l) t MWV.GKtfO^ = -(y*(<)-G; )i (c?0) ) (2-23) 

where a* = aW /(TEL 4 ) is the dimensionless kinetic factor. If we consider the 
mobilities in Scction l2.4.11 then we have the relation 



Mr 



Per 7l'(cD 



where Pc Cj i = L 2 /(TD£'), as the ratio of advection to diffusion of bulk sur- 
factants, is the bulk Peclet number and Per = L 2 /(TD-p) is the corresponding 
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interface Peclet number. The dimensionless surfactant equations with Fickian 
diffusion read as 



= 0, 



fif.c^ + cJVr.-t;.-Vr. 



Per 



Pe 



v. 



(2.24) 
(2.25) 



ai i) ^V*ci i) - I / = -( 7 :(cr)-G' M (c«)). (2.26) 



3 Phase field model 

3.1 Model for two-phase fluid flow 

In this section we will derive a phase field model for two-phase flow with surfac- 
tant generalizing the work by Abels, Garckc and Grim on phase field modelling 
of two-phase flow pQ. We start by recapitulating their essential assumptions 
and governing equations. 

For a test volume V CCI, let p denote the total mass density of the mixture 
in V and, for i — 1,2, denote by p^',Vi the bulk density and the volume occupied 
by fluid i in V, respectively. Let Ui = Vi/V denote the volume fraction occupied 
by fluid i in V. Assuming zero excess volume due to mixing, we have 

ui+u 2 = l. (3.1) 

Then the total density p can be expressed as a function of the difference in 
volume fraction ip = U2 — u±, which is a natural choice for the order parameter 
that distinguishes the two fluids, 

p (2 Hi+^) , p (2) -p (1) , p (2) +p (1) 

P = p(<f) = 2 + ~ 2 = ~ 2 V+ ~ 2 • 

As in [I] [55] , we assume that the inertia and kinetic energy due to the motion 
of the fluid relative to the gross motion is negligible. Therefore we consider the 
mixture as a single fluid with velocity v. If one chooses v to be the volume 
averaged velocity then the prototype diffuse interface model for incompressible 
two-phase flow with different densities is: 

V-t> = 0, (3.2) 
d t (pv)+V-(pv®v)=V-T, (3.3) 
dtip + V-(<pv) = -V-J v , (3.4) 

where T is a tensor yet to be specified, J v is a flux related to the mass flux J 

by 

(pW~p^)J v = 2J. (3.5) 
As a consequence of (|3.4p we obtain the mass balance law 



<9 t p + V • (pv) = -V • J. 



(3.6) 
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Our goal is now to extend this model to the case where surfactants are 
present, distinguishing the cases of dynamic and instantaneous adsorption. We 
proceed as in the sharp interface setting by postulating appropriate mass balance 
equation(s) for the surfactant and deriving models from constitutive assump- 
tions such that thermodynamic consistency is guaranteed. 

3.2 Dynamic adsorption (Model A) 
3.2.1 Mass balance equations 

We will use the distributional forms for the bulk and interfacial surfactant equa- 
tions to derive the phase field surfactant equations. Since the sharp interface 
is replaced by an interfacial layer, we consider rcgularisations of XfiW an( i <^r 
that appear in (|2 . 14(1 . (|2.15|) . (|2.16|) . In the context of phase field models, many 
regularisations of the delta function are available from the literature [3H1 [THl 03] , 
but it will turn out that the Ginzburg-Landau free energy density 



is a suitable regularisation for a multiple of 5r, where e is a measure of interfacial 
thickness and W(ip) is a potential of double- well or double-obstacle type [8] 
with equal minima at p = ±l and symmetric about p = 0. For example, one can 
choose W(ip) = \ (1 — ip 2 ) 2 for a potential of double- well type or 



for a potential of double-obstacle type. However, in the following derivation 
we assume a smooth potential for convenience. The potential term W{p) in 
S(p,\7<p) prefers the order parameter p in its minima at ±1 and the gradient 
term |Vy>| penalises large jumps in gradient. This leads to the development of 
regions where tp is close to ±1 which are separated by a narrow interfacial layer. 
For the regularisation of Xn( 2 h we consider ^{p) to be a non-negative cut-off 
function such that £2(1) = !, 1) = and £2 varies smoothly across |y|<l. 
For example, in some of the subsequent numerical experiments we used 



Similarly, £i(ip) = 1 — £,2^) will be the regularisation of Xnw ■ 

Our ansatz for the case of dynamic adsorption of the surfactant to the in- 
terface is motivated by the distributional formulation in (|2.14j) - (|2.10l) 



5(ip,Vp) = ^\Vip\ 2 + -W(p) 





dt & c« ) + V • (£ Sv ) + V • (£ J W ) = Sji , i = 1 , 2 , 
d t {K8c T ) + V • (KSc r v) + V ■ (K6Jr) = S(ji +32), 



(3.7) 
(3.8) 



where Ky^O is a calibration constant which depends on W, chosen such that 
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K6(ip,V<p) regularises Sr, see [32]. In particular we set 




W{Unh(z/V2))dz = 2^2/3, for the double-well, 
2W(sm(z))dz = 7r/2, for the double-obstacle. 



Furthermore, Jc is the bulk surfactant flux, Jr is the intcrfacial surfactant flux 
and ji,i = 1,2, denote the mass exchange between the bulk and the interfacial 
regions. In the above prototype model we allow the situation where there are 
surfactants present either in both bulk phases or in just one bulk phase. We 
denote the former as the two-sided model and the latter as the one-sided model. 
In the one-sided model, we set = 0,£i(</?) = 0, ji = 0, J c .\ = and we drop the 
subscripts so that equations (|3.7[) , (|3.8[) are written as 

d t (l<5{v,Vtp)c v ^ +V ■{K5{tp,VL P )c T v)+V ■(K5{tp,VLp)Jr)=-5{v,Vtp)j. 
Observe that, for a test volume V(t) with external normal v, we have 

4(E/ ^+1 K5A=-[ {hjp+bJ<n+K5J T ).u, 

i^i,2 Jv Wl J v(t) ' JdV(t) 

which is analogous to (|2.2[) . 
3.2.2 Energy inequality 

We introduce a Hclmholtz free energy density a((p, Vy?,cW,c r ) which will play 
the role of the bulk and intcrfacial free energy density for the diffuse interface 
model. As in the sharp interface setting and in analogy to (|2.3p the total energy 
in a test volume V is the sum of the kinetic and free energy: 

l^e(v,ip,Vip,c^,c r )^jy^f+ J a(<p,Vtp, C W,<F) (3.9) 

where 

a( V yv,c,c r ) = KS&,Vvh(c r )+Z 1 (v)Gi(c^)+Uv)G2(c {2) ). 

Since K5((p,V(p) approximates 6r we can consider the first term as an approxi- 
mation of the surface free energy density. We assume that the free energy densi- 
ties satisfy 7" > 0,G"/ > and that the following dissipation law holds pointwise 
in V: 

-V:=d t e + V-(ve) + V-J e <0 (3.10) 

where J e is an energy flux that we will determine later. 
From (|3.6p and Q3.3P we have 

d t ( £ ¥-) + V • (^") = -^f V J + (V -T) ■ v + [(V -J)v] ■ v 

= — ^V- J+(V-T)-i>+[V-(v<8> J)]-v-[(J-Vv)]-v 

= V- (-^J + T-t-vj-T: Vv + [V-(v® J)]-v 

= V- (- {T x + [v(3~J] ± )v S J - (T+ (w® J)) : Vu. 
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We use the identities 

djV(p = V<9 4 V- (Vv^Vtp, d'(ab) = ad'b+bd'a, 
d' t {K5{ Vl V v ) 1 {c T )) = d' t {K5) 1 (c T ) + 1 '{c T )d' t c T K5 

= d:{K5) 1 (c r )+ 1 l {c T )d' t {K5c T )- 1 \c T )c v d' t (K8) 1 
d' t (d {<p)Gi (c« )) = d' t (fc (^)cW )G' i (c« ) + 9 t * (6 (<p))(Gi (c« ) - c« G- (c« )) 
to obtain after some lengthy calculations that 

i=1.2 

+ V ' ( E £ (V) ( G * ( c(l) ) - G * ( c(4) ) cW ) - V • (KeaV^p) + -aW'^))) 

i=l,2 

+ KSJr ■ VV (c r ) + 6 ■ VG; (c«) + 6 J< 2) ■ VG' 2 (c< 2 > ) 
- 5ii( 7 '(c r ) - GU C (1) )) - 5j 2 (7'(c r ) - G^)) 

+ J V -V( £ e, , M(G l ( C (i;) )-G:(c«)cW)-V.(^ £( 7V^) + - ( T^'(^)) 

i=l,2 

+ (V -v) (- y( Yl & ^) ( G » ( c( '° ) ~ G * ( c<4) ) cW ) - V ' (^V V ) + f <rl*») ) 

i=l,2 

+ (V • v) (K8o + UGi(c^)-G' 1 (c^)c^)+UG2 {c {2) ) - G 2 (cW ) C W )) 

-Vv: (T + txg> J + A'ecrV^V^). 

In the case where the surfactant is present in only one of the bulk phases, a 
similar calculation shows that we obtain the above form for —T> without any 
terms involving the subscript 1. 

In any case, we choose J e so that the divergence term cancels. 

3.2.3 Constitutive assumptions 

We set 

H = -V • {Kea(c v )V V ) + ^a{c T )W^) + £ £^)( G i(c (i) ) - G'MV } ) 

t=l,2 

and make the following constitutive assumptions: 

Jr = -M r (c r )V 7 '(c r ), 

J 2 =/3 W (7'(c r )-G:( C «)), 
J v = -m(ip)X7fj, 

for some non-negative function m(ip), and the /?W,i = l,2 are given by 
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We choose the tensor T to be 



T 

i=l,2 

- v ® J - A'ecrV^® V<y? + 2r)(tp)D(v) -pi 

where p denotes the unknown pressure, 7](tp) > denotes the viscosity and from 
(|3.5|) the volume diffuse flux J is given by 

J=- P (2) -P (1) m ( y )VM. 

Since the interface thickness will be of order e it turns out that the term 

V • (Ka(S(ip, Vtp)I-eVtp® Vip)) 

scales with e~ 2 , while the term 

V.(a(G 1 ( C ( 1 ))-G' 1 (c( 1 ))c( 1 ))J + 6(C? 2 (c( 2 ))-^( C ( 2 )) C ( 2 ))/- W /) 

scales with e~ l , the same order as the pressure p. Hence we absorb the latter 
term as part of the pressure and reuse the variable p as the rescaled pressure, 
leading to 

T = Ka(c r )(5(tpytp)I-eVtp®V(p)-pI + 2Ti(tp)D(v) + v® J. 

We remark that the term V- (Ka5((p,^(p)I) in the momentum equation is re- 
quired to recover the surface gradient of the surface tension in the asymptotic 
analysis. It is present also in other diffuse interface models with Marangoni 
effects [321 ED III] • 

With the above assumptions we obtain the energy inequality 

-V = -m(<p)\Vfi\ 2 - M«(cW)6(^)|v^(c«)| 2 -2r ? ( ¥ >)| J D(t;)| 2 

i=l,2 

-/3 W (5(^,V^)|7'(c r )-GKc W )| 2 - KM r {c r )8{ip 1 V V )\V 1 ' {c r )\ 2 <0, 

and the diffuse interface model (denoted Model A) for the case of dynamic 
adsorption reads 

V-« = 0, (3.11) 

dt(pv) + V- {pv® v) = V • ( -pI+2r](<p)D(v) + v® P W ^ m((p)V fi\ (3.12) 

+ V ■(Ka(c r )(S((pyip)I -eV(p®Vip)), 
<9 t V = V-(m(<p)V/Li), (3.13) 

M + V ■(Kea(c T )V<p) = -a(c T )W'(v) + £ ^{G^c®) -^(c«)c«), 

£ z — ' 

i=l,2 

(3.14) 

d' t &)c {l) ) = V • (M« (c« MtfVG'i (c« ) ) (3.15) 
+ /3«^,V^)( 7 '(c r )-GUcW)), i=l,2, 
d*(KS(ip,Vip)c r ) = V ■ (^M r (c r )KS(ip,yip)V-f'(c r )j (3.16) 

-%,V^)^/3«( 7 '(c r )-G:( C «)). 
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3.3 Instantaneous adsorption, one-sided (Model B) 

To model instantaneous adsorption, we assume that the surfactant is insoluble 
in one phase Similar as in Section l2.4.2l we assume that the bulk surfactant 
in fj( 2 ) and the interface surfactant are in local thermodynamical equilibrium. 
This means that the bulk chemical potential G'(c?>) and the interface chemical 
potential 7'(c r ) are equal. Hence we impose the constraint 

7 '(c r ) = ^(c( 2 )) 

in order to replace c r . For this purpose, since 7' is strictly monotone (recall 
that 7 is strictly convex) we may set 

9 (cW) = ( 1 r 1 (G> 2 (cV)) = C r. 

We then consider one surfactant mass balance equation which we obtain by 
adding (J37f]| for i = 2, d^SJ) and setting jj = 

5 t *(^)c+A'%,Vv)5(c)) + V-(f( V )J c +A'%,Vv)Jr) = 0, (3.17) 

in place of p.7[) and f|3 . 8[) (for convenience, we drop the index 2 of £2, c^ 2 \ Jc 2 "* 
etc.). 

The energy of the system is given by 

e{v,i P ,Vi P ,c) = -p\v\ l + K8{ip,Vip) 1 {g{c))+i{ip)G{c), 

and we set 

M = -V • (Kea(g(c))V<p) + * a (g(c))W (G(c) - G'(c)c), 

where 

a(g(c))= 1 (g(c))~ 1 i(g(c))g(c)= 1 (g(c))-G , (c)g(c). 
Then, a similar computation as in the previous model yields the following 

-V = V- ( J e -J^f + {v®J)v- K5j'(g(c))J T - ^G'(c)J c + Kea(g{c))Vipd'ip) 
+ V- (t ± v + J v ^ + J v -Vn-Vv: (T+v®J + Kea(g(c))V(p®'V<p) 
+ K5J T -Vi{g{c)) + iJ c -VG'{c). 

We choose J e ,T,J v as in Model A. Furthermore, we assume that 
J c = -M(c)VG'(c), J r = -M r (g(c))Vi'(g(c)) = -M r ( 5 (c))VG'(c). 
We then get the energy inequality 

-V=-2 v &)\D(v)\ 2 - m &)\V f i\ 2 -(M(c)av)+K6(vy^M r (g(c)))\VG / (c)\ 2 
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The diffuse interface model for this case (denoted Model B) is 

V-v = 0, (3.18) 

d t {pv) + V -{pv®v)=V ■ (^-pI + 2i]((p)D(v)+v® 77(21 m(v?)V/ij 

(3.19) 

+ V ■ (Ka(g(c))(6(cp,Vcp)I -£Vcp®Vcp)), 
%<p = V-(m{<p)Vn), (3.20) 

A t + V • (if £ a( ff (c))V^) = j*(g(c))W'(<p)+?(<p)(G(c) - G'(c)c), (3.21) 

dt(Z(<p)c + K8(<p,V<p)g(c))=V-(M(c)t(<p)VG f (c)) (3.22) 

+ V-(M r ( 3 (c))W(^V V )VG'( C )). 

3.4 Instantaneous adsorption, two-sided (Model C) 

We now derive an alternative model for instantaneous adsorption that is two- 
sided. Since we assume local thermodynamical equilibrium, the chemical poten- 
tials Gi(c^) ,G' 2 (c^) and 7'(c r ) are equal on the interface. We hence introduce 
a chemical potential, denoted by q, and consider this as unknown field rather 
than the densities of the surfactants. Since the free energies Gi,7 are strictly 
convex, their derivatives are strictly monotone and we obtain a one-to-one cor- 
respondence between the and q, i.e. 

C ( 1 ' = (G' 1 )- 1 (g), C W = {G' 2 )-\q), c r = ( 7 ')- 1 (?). 

We then also may write the surface tension as a function of q, 

a(q) = <T(c r (q))=l(c r (q))~c r (q)q. 

Summing (|3.7[) for i = l,2 and (|3.8D we obtain the conservation of surfactants as 
follows: 

d:(^)c^( q )+U^c {2) (q)+K6^y^c r ( q )) 

= -V • (&( V )J« +Uv>) 4 2) +K6{<p,V<p)J r ). 
The energy density of the system is given by 

e(^,Vv 3 ,^, g ) = ip|^|V6MGi(c( 1 )(g))+6MG 2 ( C ( 2 H9))+^,V^)7(c r ( g )) 
and similar computations as in the previous models yield 

-V = V • ( J e - J lj f + (v®J)v - KdqJr -tiqjP- &qjP + Kea(q)Vtpd;<p) 
+ V-(T- L v + J v fij +J v -Vfi-Vv: {T + v®~J + Kea{q)V(p®V(f>) 
+ KSJ T • Vg+6 {ip) J c (1) • V?+6 0) J c (2) • V?, 
where 

M= E ^(v)(Gi(cW(?))-?cW(g))-V-(lir e cr(g)V V ) + -cr( g )W( V ). 
z — ' e 

i=l,2 
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Choosing J e ,T,J v as before (but with the c^' now as functions of q), and 
setting 

J c « = -Mp( C W(?))V ?) Jr = -A/r(c r ( (? ))Vg, 
leads to the following energy inequality: 

-P = -2^)|^)| 2 -rnM|V M | 2 

]T M^(c^(q)M^ + M r (c r (q))K5(<py<p)) \Vq\ 2 <0. 



The diffuse interface model for this case of instantaneous adsorption based 
on the chemical potential as a field (denoted Model C) is 

V-w = 0, (3.23) 
d t {pv)+V- {pv®v +pl - 2r}(<p)D(v)) = V • ( - v ® p (2) -p (1) m(y)V/x) (3.24) 



2 

+ V- (i^5-(g)(<5J-eV<p®V(/5)), 
a t V = V-(m(^)V/i), (3.25) 

M + V-(^( g )V^)-|a(g)W^)= 2 eKv)(G(c«)-gc«), (3.26) 



i=l,2 



a;(^ 1 \q)+^ 2 \q) + K5 < F(q)) = ]T V-(Af c W( C W(?))6V 9 ) (3.27) 

i=l,2 

+ V-(A/r(c r (g))/«Vg). 

3.5 Specific models 

3.5.1 Insoluble surfactants 

Similar as in Section 12.4.31 we can consider a phase held model for insoluble 
surfactants. The resulting model is a system for the unknowns v,p,(p,p,,c T and 
is obtained by setting & = and = in ([3"TT]) - (|3TTo) . 

3.5.2 One-sided model with non-instantaneous adsorption 

It is also possible to consider a one-sided version of Model A by setting £i = 
and neglecting the unknown 

3.5.3 Mobility for the phase field equation 

We will choose the functional form of the mobility to be 

m(ip)=m 1 (l-ip 2 ) + , 

where mi >0 is a constant and (■)+ denotes the positive part of the quantity 
in the brackets. This degenerate mobility switches off diffusion in the bulk 
phases away from the interfacial layer. In this case, the phase field equations 
(|3.13[) , (|3.14p lead to a pure advection of the interface, see [T]. 
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3.5.4 Diffusivities 

If we set 

M {l) =D (i) Mr(c r ) = Dr 

c c G'l{c^y r[ ' r 7 "(c r )' 
for constants and Dr, then we derive Fick's law for the surfactant 
J® = -j)(')Vc(' ;) , J r = -L> r Vc r . 

3.5.5 Obstacle potential 

If W is chosen to be a potential of double-obstacle type, then equation (|3.14p is 
formulated as the following variational inequality: For all ifj £ K, := {rj £ H 1 (fl) : 

W<i}, 

K 

-p{ip - + Aeo-(c r ) V<p • (V^ - Vy>) H CT(c r )iy' Up) (-0 - 

e 

^^^)(Gi(c«)-GKcW)cW)(^-^)>0. (3.28) 



3.5.6 Reformulation of the momentum equation 

A short computation shows that 

pVtp = V • (Ka(5(<p, Vtp)I - eVtp <g> V</?)) - if V<p) Ver 

i=l,2 

hence the momentum equation (|3.12p can be reformulated as 

d t {pv) + V -(pv®v) = V- ( -pJ+2r?((^)£)(t;)+i;(8) p (2) -p (1) m(y)V^j 

+ /jV<^ + A'<5(^,V< / ?)Vct- ^ ^(93)(Gi(cW)-G;(cW)c«)V93. 

i=l,2 

3.5.7 Non-dimensional evolution equations 

We consider the following dimcnsionless variables 

e m(<g)E \iL 

d* = Ld, = —, m* = , /j* = — 

with the characteristic length L, the scale £ for the surface tension and a char- 
acteristic velocity V. In addition we scale the bulk densities by C, the interfacial 
density by C r and similar to the density, the viscosity T)Up) can be decomposed 
to ri = uiri^ + v,2T)^ 2 ' ■ The dimensionless density and viscosity are 

p* = p/p {2) =u 1 X p + u 2 , 77* =r)/rf 2 ' =ui\ v +u 2 
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where A p = p^ / p^ , A,, = jrf 2 ^ are the density and viscosity ratios. Set Re = 
{-p {2) L 2 )/{T^), Ca= (?/ 2 ).L)/(TE) to be the Reynolds and capillary numbers 
respectively. Then the dimensionless fluid and phase field equations are 

V* •«* = (), (3.29) 

V Re 2 

-V* • (lf<7*(<5*J — e^V*^® V»v?)) , 



ReCa 

(3.30) 

c\" <^ = V* • (m*(< / 9)V*/^*), (3.31) 
V*-(iY£^V^)-- ( 7^ , (^) = - Ai>! ,+ V ^'(^(G^.z^Vg'.^cI )^), 

(3.32) 

where p 1f = {pT 2 )/{L 2 p^ 2 ^) is the rescaled pressure. The reformulated momen- 
tum equation from Section 13.5.61 has the dimensionless form 

/ 277* 1 — A p \ 

dt, (p*v*) + V* ■ (p*^* ® f *) = V* • + — -D(t>*) + ® — - — m*(( ( 9)V*^*J 



+ 7rirfM*V^ + /a„V* ( T st N ) (3.33) 
ReCa V / 



i=l,2 

The dimensionless surfactant equations for Model A arc 

^.te^-V*- (mW&V.G^c? 5 )) =^ i) <y.('/,(cJ)-G^(cl i) )) J (3.34) 

(jf^cf)-v..(iirMr,.^v*7i(^)) = -<y* y /^WC^-cgCc^)), 

i=l,2 

(3.35) 

where = /3ETL 4 . For Model B, the dimensionless surfactant equation reads 



9 t * (£c*+K5*g*) - V* • fM c ,4V*G;(c*)+irMr,*5*V*G;(c*)J =0, (3.36) 
and for Model C, it reads as 
^(6cl%0+6cl%*)+A^(9*)) 

-V* • (Af$£iV,g. +A/i, 2 . ) 6V,g, + iCMr,»i»V»g») =0. (3.37) 



If we consider the mobilities in Section [3. 5. 41 the dimensionless surfactant equa- 
tions for Model A arc 

^(fi^Vv.^^^c^^^^Cyj^)-^^)), (3.38) 

a£ -v.-(^<y.v,£) y -g^)). (3.39) 

r i=l,2 

For Model B, the dimensionless surfactant equation with Fickian diffusion reads 
<9 t * (£c+jr£*5*)-V*-(^fV*c* + — J*V*c*) =0. (3.40) 
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4 Sharp interface asymptotics 

In this section we identify the sharp interface limit of the diffuse interface models 
introduced in the previous section by the method of matching formal asymptotic 
expansions. The procedure is based on the assumption that there exist a family 
of solutions, sufficiently smooth and indexed by e, to the diffuse interface models. 
For small e, we assume that the domain 57 can at each time t be divided into two 
open subdomains f2 ± (t;£), separated by an interface T(t;e). Furthermore, we 
assume that the solutions have an asymptotic expansion in e in the bulk regions 
(away from T(t;e)) and another expansion in the interfacial regions (close to 
r(t;e)). The idea is to analyse these expansions in a suitable region where they 
should match up. We will apply this method to Model A, where we distinguish 
two different scalings of aW, namely 0(1) and 0{e). In the last section we 
briefly outline the procedure for Models B and C. Details of the method can 
be found in [551 UJ f° r the smooth double-well potential and in [SJ [7] for 
the double-obstacle potential. We remark that for some specific models this 
procedure has been rigorously justified (see [2l IT5l ITT] ) . 

4.1 Outer expansions, equations and solutions 

We assume there exist the following asymptotic expansions in e for u e = 
u(t,x;e) G {v e ,p e ,tp e ,fj, e ,Ce ,c£} in the bulk regions away from the interface 

u s (t,x)=u(t,x;e)=u {t,x)+eu 1 (t,x) + O(s 2 ). (4.1) 

Substituting these expansions into Model A and (|3.14j) to order —1 gives 

Q = a{cl)W'{ip ). 

As cr>0, we obtain the stable solutions ipQ = ±l. We denote and f^ 1 - 1 to 
be the sets where tpo = 1 and ipo = — 1 respectively. 

The zeroth order expansions of the fluid equations yield 

V-«o = 0, 

dt{p (i) v ) + V • (p^v ® v - 2r)^D(v ) +p I) = 0. 

The bulk surfactant equation gives to the zeroth order 

d t cf +v ■ Vc« - V • {Miicf^G'^)) = 0, i = 1,2. 

Observe that 6(ipo, V<^o) = so that (|3.16p fully degenerates in both domains 
f^ 2 -* and whence Cq remains undetermined in the bulk. Similarly, /iq is 

undetermined in the bulk due to the degenerate nature of the mobility m((po). 

For the double-obstacle potential, equation (|3. 14|) is replaced by (|3 . 28[) which, 
to order —1, is the variational inequality 

[ a($)W'(cp )(iP -ip )>0, VV^oe/C. 
Jn 

Here, W'(tp)=—tp + dL_i i]((p) where dl is the sub-differential of /r_i ii. Then 
the above can be expressed as 

cr(Cg )</?oC0o - tp ) >0, V^qG/C. 
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Since a > 0, this implies that ipo must take the values ±1 and we can define sets 
as in the case with the double-well potential. 

4.2 Inner expansions and matching conditions 

Let us assume that the zero level sets of tp s converge to some hypersurface T 
moving with a normal velocity denoted by ur as e— >0. Close to T, we denote 
by d(t,x) the signed distance function of a point x € f2 to T with the convention 
d(t,x)>0 if x£Q( 2 '(t), and set z{t 1 x) = d{t,x)/e. We write each field u(t,x) 
close to r in new coordinates U(t,s,z) where s are tangential spatial coordinates 
on r. The upshot is 

d t u = -\u v d z U + dtU+ h.o.t., 
V x u=-d z Uv + V T U+ h.o.t., 

£ 

A x u= \d zz U --k8 z U - z\S\ 2 d z U + A r U + h.o.t., 

e 2 - e 

where v = \7 x d is the unit normal pointing into fi' 2 ), d^(-)~dt(-)+ur^-V x (-) 
is the normal time derivative, Vr is the spatial surface gradient on T, n is the 
mean curvature, |<S| is the spectral norm of the Weingarten map S, Ar is the 
Laplace-Beltrami operator on T and h.o.t. denotes higher order terms (see the 
appendix of [T] for a proof). 

We assume that the inner expansions of unknown fields u G 
{v e ,p £ ,ip e ,fi e ,ci l \c^} take the form 

u(t,x;e) = U(t,s,z;e) = U a (t,s,z)+eU 1 (t,s,z)+0(e 2 ) 

with inner variables Ue {V,P,$,M,CW,C* r }. We assume that $ satisfies 

$(t,s,0;e) = 0. 

Regarding the double-obstacle potential, we further assume that $ is monotone 
increasing with z and the interfacial layer has finite thickness of 21, where the 
value of I will come out of the asymptotic analysis (see [5]). For the double- well 
potential we take l = oo. Furthermore, we assume that 

$(t,s,l;e) = l, $(t,s,-l;e) = -l. (4.2) 

In order to match the inner expansions valid in the interfacial layers to outer 
expansions we employ following matching conditions [2S]: As z— >±Z, 

U (t,s,z)~v£(t,x), (4.3) 

d t U (t,8,z)~Q, (4.4) 

0,l7i(f,*,*)~Vuj} : (f,!B)-i/, (4.5) 

d z U 2 (t,s,z)~Vuf(t,x)-v+((vV)(vV)ut(t,x))z, (4.6) 

where denotes the limit lim^o uq (x ± 5v) at a point a;er. 

If the bulk fields are not determined by any equation, i.e., if u = c r or u = fj,, 
then we assume that the derivatives of the inner expansion in z remain bounded 
as z— >±l. More precisely, we assume that 

d z Co, d z C\, d z C2, d z M , d z Mi are bounded as z->±I. 
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Moreover, we assume 

SzCg 1 ^, d z Ci are bounded as z — y+l, d z C^\ d z C[ 2 ^ are bounded as —I, 

since c' 1 - 1 is not defined in fl^ and cJ 2 ^ is not defined in fl^K Similar assump- 
tions arc made for the asymptotic analysis of Models B and C. 

4.3 Asymptotics for Model A 

We begin by stating a few expansions of the most complicated terms for later 
use. These can be obtained by some short calculations. First, 

eV-(<7(c r )Vip®V v >) = \d z (a(c r )(d z $) 2 v) + -d z (a(c r )d z <f>Vr<S>) 

+ iv r -(cr(c r )(^$) 2 iv(g)i/) + Vr-(cr(c r )92$(iv(X)Vr$ + Vr$«)^))+ h.o.t. 

where Vr - of a 2-tensor is the surface divergence applied to each row. Then, 
setting £(A) = \{A + A ± ) for a tensor A one can show that 

V • (v(<P)D{v)) = \d z (r)($)£(d z V®»)u) 

£ z 

+ -d z ( V ($)£ (V r V» + -Vr ■ (d z V® v)) + h.o.t. 

e e 

Finally, observe that 

8(tp,Vtp) = ^ |d 2 $| 2 + - + 1 |V r $| 2 + h.o.t. 

Z£ £ 2 

4.3.1 Inner equations and solutions to leading order 

The order —3 terms in Q3.16P give 

K 3 Z ( M r ( Cl ) ( 1 1 8 Z <£ o | 2 + W ( $ ) ) d z 7 ' ( C r ) ) = . 

Integrating from -i to z and matching conditions (|4.3[) and (|4.4|) applied to $ 
yields 

A/r(C r )(i|9 z <I>o| 2 + VF($oW))O z7 '(C7 r (z)) = 0. 
We conclude that 

d z i{Cl{z))=Q whenever |$ | < 1. 

Since 7" > 0, we obtain that 

d z Cl (z) = whenever |<fr | < 1 

which means that Cq is constant across the interfacial layer. Since the surface 
tension is given by ct(Cq) = j(C^ ) — Cg j'{Cq ), we also obtain 

d z a(Cl (z)) = whenever |<J> | < 1- 
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To order —1 in (|3.14[) we have 

AV(C r )(-S zz $ o + H"($„)) = 0. 

We can choose <I>o such that it is independent of s and solves 

-0„*o + W(*o)=O, (4.7) 

with $ (0) =0 and $ o (±0 = ±1 - with thc double-well potential W(tp) = |(1- 
ip 2 ) 2 we have the unique solution 

<J>o(z) = tanh(z/\/2), 

while for thc double-obstacle potential, a unique solution to 

-a 2 ,$o-$ = 0, |$o|<l, $o(M,0)=0 

is 

{+1, forz>§, 
sin(z), for |z|<§ , 
-1, forz<— f, 

so that I = \ and from (J472J) we deduce that 

$i(t,«,±f) = 0. (4.8) 

Multiplying (|4.7[) by 9 2 $0j integrating from — / to z and applying matching to 
$o yield the cquipartition of energy 

\\d z ^(zt = W(^(z)). (4.9) 

The order —1 term in the mass balance (|3.11[) gives 

(d z Vo)-u = d z (V -v>) = 0. (4.10) 

Integrating from — I to I and matching (|4.3[) applied to Vq imply that Vq-v is 
constant in z and 

v^-u~ lim Vb(z)-i / = lira Vo(z) ■ f = i'q 1 '' • u, (4-H) 

i.e., the normal velocity is continuous across the interface. 
Equation (|3 . 1 5[) gives to order —2 

d z (M« (C« ($o)Gr (cf )d z C® ) = 0. 

In the two-sided model, for i = 2 we integrate from — I to z to obtain 

M 2 (cf )6(* (s))G£(cf (z))d z cf (z) = 

as ^2( — 1)=0. Since G'{ >0 we have that 8 z Cq =0. Similarly for Cq where 
we integrate from 2 to +/ to obtain 

M 1 (C^)^ o (z))G r {(C^\z))d s C^ ) (z) = 
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as £i(+l) = 0. Thus 8 z Cq 1 — follows from the same argument. In the case of 
the one-sided model, we argue as above to obtain d z C'o = 0. 
Equation (|3.13[) gives to order —2 

O = z (mi(l-*jj)+0,M o ). 
Integrating from — I to z and matching (|4.3|) applied to $o gives 

= mi(l-$g(«)) + 9,M («). 

For |$o| < 1 we have d z Mo = 0, hence the term V- (u<g> ^ ( } ) m(y)V/x) plays 
no part in the order —2 expansion of the momentum equation (|3.12[) . To leading 
order the momentum equation gives 

= 2dM*o)d, V ). (4.12) 

With the usual trick of integrating with respect to z from — I to a limit denoted 
by z again and applying (|4.4[) to Vq we obtain that ry($o)^ Vo = 0. Since 77 > we 
conclude that d z Vo = so that, using (|4.3|) . the tangential velocity is continuous 
across the interface: 

[«o]? = 0. 

4.3.2 Inner equations and solutions to first order 

Equation (|3.11[) of the mass balance yields to zeroth order 

d z V 1 -v + Vr-V o = 0, (4.13) 

while equation (|3.13[) gives to order —1 

(-«r + V o -i/)0,$o = 0,(mi(l-*g)+a,Mi), 

where we used that d z M$ =0. Integrating from —I to +1 and applying (|4.3[) to 
i»o and (|4.5[l to Mq (see also the remark after further down after (|4.2p ) then 
imply that 

2(u r - t>o • v) = [mi(l - #g) + a,Mi]l{ = 0, 

and we obtain 

u r = v -v. (4-14) 

Using equipartition of energy (|4.9[> . B z Cq — and «r = t>o • we obtain from 
(|3T5)l at order -1 

2 J 9W(y(C'f)-GU^))W($o) = -5.(M i (C<^)C i (*o)^(Gr(^)C'«)). 

(4.15) 

In the two-sided model, for i = 2, integrating (|4.15|) from —I to +1 and using 
P~5j) leads to 

= [M 2 (cf )6($o)G' 2 '(Cf )fl a cf ]t{ + 1 + 2/3( 2 ) ( 7 '(C r )-G 2 (G^ ))VF($ )^ 
= A/( 2 )( c ( 2 ))VG 2 (c^)^ + /3 (2) (7'(^)-G 2 (c^))W. 
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Proceeding similarly for i = 1 and we recover the following free boundary con- 
ditions 



The argument for the one-sided model is similar to the above case with i = 2. 

Using 8 z Cq = 0, ur — vq ■ v, and the equipartition of energy, after integrating 
from —I to z and matching, equation (|3.16[) gives to order —2 

2M T {Cl)W(^{z)) 1 "{Cl{z))d z C{{z) = Q. 

Since 7" > we have that 

d z C\ = whenever |<j> | < 1. 
Equation (|3.14[) for the chemical potential gives to zeroth order 
M = Ka(C$){-d„$i + ($ )*i) +Ka'{Cl)C\ (-d 22 $ + W'(* )) 

S * ' 

=0 

-A"V r -(a(C7 r H5 2 <l>o+ £ ^(G^C^-G^C^C^). 

i=l,2 

To obtain a solution $1, a solvability condition has to hold. Multiply the above 
by dz&o and integrate from — I to +1 gives 

/•+2 r+l 

J ^ M d z <P dz = K J ^ <r(C£)(-0„*i0,# o + W ,/ (#o)*i0.*o)«fe 

-K [ +l V r -(a(C^)(d z <P ) 2 dz+ [ + J2 a^){G t {C^)~G' t {C^)C^)d z ^ dz. 
J —i J-i ■, n 



Integrating by parts, using <9 2 Cq =0 7 d z C^ = and matching lead to 



2/x = K f a(C r ) (d zz <S> ~ W'(<$> Q ))d z 3> 1 dz~K[a{Cl){d z <$> Q d z <$> 1 - W"(* )*i)] 
" , ' 

=0 

-KVr-{o{C T )v) [ + (d z <Z> ) 2 dz + ]T (C^ ) - G\ (C^)C^ )& ($0)] _j ■ 



We use the fact that W'(±l) = for the double- well potential and (|4.4p to cancel 
the first jump term. Furthermore 

V r • (<r(Cg» = a{Cl)V T ■ v + V r <7(C r ) ■ v = -k<t(C£) 

v ' 

=0 

and by equipartition of energy f|4. 9[) we deduce that the solvability condition is 
2^a{cl)n+[G l {c i i ) )-G' l {c ( i ) )c ( i ) ]l (4.17) 
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For the double-obstacle potential, the equation for $1 is expressed as a varia- 
tional inequality: For all ipo€zfc 

A'(-a(C r )(a M $ 1 +$ 1 )-a'(C r )Cr(^$ + $o)-5 2 $oVr-(a(C r )^),V'o-$o) 
> (m - e-($o)(G l (^ ) )-G' J (^ ) )C^),V'o-<l>o). 

»=1,2 

Whenever |$o|<l, testing with ipQ — $o + ipo with either a non-positive or a 
non-negative ipo we obtain the equality 

- M - AV(C£) (5 22 $! + <Di)-&' (Cg")Cf (d 22 $ + *i) 
-A^ 2 $ Vr-(a(C >)+ £ ^(*o)(a(C«)-^(C»)C«) = 0. 

i=l,2 

Multiplying by 9 2 $o an d integrate from — I to gives after matching 
2 / i -a(c >- ]T fe(¥'o)(^(c«)-G ! K4 <) )c«))]l| 

i=l,2 

=/v / - ( x(c r )(a 22 $ 1 +$ 1 )a 2 $o^ 

r+l 

= -A[a(C r )(9 2 $ a 2 $i+$ $i)]lS+AV(C r ) y ^ 5,*i(fl„$o + *o)«fa. 

The last integral term is zero due to (|4.7[) . and using (|4.4D for $o and (|4.8[) 
for $i at 2 = ±Z the jump term is also zero. This leads to the same solvability 
condition as in (|4.17j) . 

Using d z M = 0, ur = v -v> 7 Vr'&o^O and cquipartition of energy, the mo- 
mentum equation (|3.12[) gives to order —1 

9 2 P ^ + ^($o)f(a 2 yi®iv) I , + r ? ($ )£(VrVo)^-c) 2 (yo® 1 ^^m($o)a 2 M lI /) 

= K |a 2 $ | 2 (V r • WC r )/)-V r • (tr(Cg>8> u)) 

where we used that Vb is constant in z. Matching (|4.5|) requires that 
]im z ^±td z V±(z) = Vd^v and hence 

d z Vi<gn> + VrVo^Vv for z->±l. 

Furthermore, a short calculation shows that 

Vr-((T(C^)I)-Vr-(<y(C^)^®u) = Vr<j(C^) + na(C^)iy. 

So upon integrating from —I to +1, matching and using that m(±l) = we 
obtain 

lpo]lu-2ri[D(v )]lv = K*(CS)v + V T *(C$). (4.18) 
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4.3.3 Inner equations and solutions to second order 

Using Ur=Vo-f, 8 z Cq = d z C\ =0 and equipartition of energy (|4.9p . equation 
(|3.16|) gives to order —1 



x(d° (2C^W^ ))+V -V r (2C%W($ )) + (Vi-v)8 z (2C r TU($ )) 

= Kd z (2M T (Cl)W{W{Cl)d z Cl) +KWr ■ (2A/ r (C r )IU($ )Vr7'(Co)) 



i=l,2 

Integrating from — I to +/, we obtain 

+i r +i 



d° t Cl + V -VvCl)K j 2W($ )dz+K J {V 1 -v)d z {2W{^ )Cl)dz 

= ^[2Mr(C r )^($o)7''(Co)<9^ 2 r ]^+Vr-(Mr(C r )Vr7'(C r ))K J ^ 2W($ )dz 
-^^(7'(C r )-G;(C«))W. 

i=l,2 

Applying the matching conditions (|4.3j) to $ an d Cq and (|4.6[) to C\ we see 
that 

[M r (C r ) (W{W{Cl)d z Cl)] + _\ = 0. 
By (HTTB")) we have that 

^ (Vi-u)d z (2C^W(^o))dz=[2(V 1 -u)C^W(^ )}t l l - J + 2d z {V 1 -v)W(<$> Q )Cldz 

r+l 

J ^ (Vr-V )2W(*o)<te = >VcJVr-t;o, 



= + c£ 



and by (j4~T^|) 

E^ W (7'(^)-GKc«))W = [jS]^. 

i=l,2 

Using 9'(-)=9 t °(-)+t)-Vr(-), we finally obtain the desired surface surfactant 
equation 

^ + ^Vr^o-Vr-(Mr( c r)Vr7'(^)) = [J c ( ^^. (4.19) 

4.4 Alternative asymptotic limit for Model A 

Let us now assume that (3^ scales with e^ 1 , i.e. 

H We 

Then we obtain instantaneous adsorption f)2 . 13|) instead of (|2.12[1 in the limit 
e— >0, which will be demonstrated in what follows. 



31 



4.4.1 Inner equations and solutions to leading and first order 

We recover [vq-v]\=Q and obtain d z V\-v-\- Vr ■ Vo = from equation (|3.11[) to 
order —1 and to zeroth order respectively. From equation (|3.13p we obtain 
O z Mq = and ur = «o " v to order —2 and to order —1 respectively. To order —2 
equation (|3.12[) gives [«o]i = 0. 

To order —3, the interfacial surfactant equation (|3.16[) gives d z C\ = 0. This 
leads to the profile $o and equipartition of energy (|4. 9[) from (|3.14p . Further- 
more, we obtain the solvability condition (|4.17[) from (|3.14p at zeroth order and 
the jump in the stress tensor (|4. 18[) from (|3.12p at order —1. 

To order —2 we obtain from (|3.15|) and (|3.16|) 

-Wd z (M«(C«)^<I> ^^ (4.20) 
Wd z (M T {Cl)K2W{^)d z {i\Cl)C T 1 ))=Y J 2VF($ )(7'(C r )-GUC«)). 



(4.21) 



Now, multiplying g^OJ) by G'^C^), i= 1,2, and (|4~2Tj) by Y(C£) and subtract- 
ing gives 

-W J2 d z {M^^ )d z G' i (C^))G' i (C^) + 2W^ ) [I'iCD-GriC^) 

i=l,2 i=l,2 

-W9 2 (A/r2A'VF($o)^(7"(C , o)Cr))7'(Co r ) = 0. 

Integrating from — I to +1, integrating by parts and using that d z Cy = yields 
= W E M c ] ti (*o) I d z G' t (C« ) 1 2 dz - W [M«& (^)d z G' t (C«)G< (Cf )] +; 

1=1,2*' ~ l 

- [M r 2KW^ ) 7 "(C^d z CW(C^t l l + £ r 2W($o)| 7 '(C r )-G;(C«)| 2 . 

i=l,2*' _i 

The third term vanishes as W(±1) = 0, and when applying (|4.4p to Gq 1 ' then 
the second term is zero, too. Hence we have 



WW Af c %(<i>o)|a 2 ^(C«)| 2 +E f 2VF($ )| 7 '(C r )-G' i (C«) 

i=l,2 1 i=l,2 1 

As all the terms are non-negative, this implies that 

d z C { o l) =0 and 1 '(C^ = G' l (C { l) ). 

4.4.2 Inner equations and solutions to second order 

Adding the surfactant equations (|3.15p and (|3.16[) . the order —1 terms yield 
2 KW($o) (9 t °C r + Vo • V r G r ) + V 1 ■ ud z {2KW(^)Cl) 

= d z (M T 2KW(^ )j'\C^)d z C^ + MrK(d z ^ d z ^ 1 + W'^ )^ 1 )j''(C^)d z C^) 
+ V r ■ (A/r2A^($ )V r7 '(G r ))+ £ d z (M^^ )G'f(C^)d z C^). 



0. 



i=l,2 



Integrating from —I to +Z and matching (|4.5p applied to leads to (|4.19p 

again. 
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4.5 Asymptotic analysis for Models B and C 

The asymptotic analysis for Models B and C are similar, hence wc will only 
sketch the analysis for Model C. In the following, the analysis for Model B can be 
recovered by setting variables with index 1 to zero and replacing c^ 2 \q),c T {q),q 
with c,g(c),VG'(c). 

First we express (|3.27p as 

^(a(^c (1) (g)+6Mc (2) (<?) + ^( V3 ,V^) C r ( (Z )) + V-J = 0, (4.22) 

where 

J:=-(MP^) + M^^) + MrKS(ipy^))Vq. 

Based on the outer and inner expansions of 5((p,V(p), we assume that J has the 
following outer and inner expansions: 

J = e- 2 J^ lk + e- 1 J b f k + J bulk + . . . , 
J = e~ 2 J™1 + e- 1 ^ + Jj nt + . . . , 

where, for example, 

jbuik = Q) jbuik = -M T {(F )W{ VQ )Vq Q , 
J™ = -M r {Cl){\ |9 2 $ | 2 + W{$ ))d z Q Q u. 

The matching conditions for J are as follows (see [H]): As z^±l, 

J^(t,s,z)~0, d,J<%(t,s,z)~0, (4.23) 
Ji!\ t (t, S ,z)~(J bl ; lk ) ± (t, a; )- I v, d z J™l(t,s,z)~0, (4.24) 
J^(t,a,z)~(J^) ± (t,x) + V(J^ k ) ± (t,x)'VZ. (4.25) 

4.5.1 Outer equations and solutions 

From equation (|3.26|l we obtain to order —1 

= a(q o )W'(<p o ), 

from which we obtain stable solutions ipo = ±1 and regions H,^,^ 2 ^ defined as 
in previous models. We also recover the usual fluid equation, incompressibility 
condition to zeroth order. 

With respect to the surfactant, to order —1 we have 

J^ = -M T KW{ip yjq = Q. (4.26) 

To zeroth order we recover the bulk surfactant equations from (|4.22[) : 

a;(a(^o)c (1 H<?o)+6(^o)c (2 ^go))-V-(A/i 1 )a(^o)Vg + M c (2) 6(v'o)V (Z o)=0 

where £ 1 (( Po ) = £ 1 (-l) = Q in fi< 2 ) and £ 2 (^ ) =6(1) = in 



33 



4.5.2 Inner equations and solutions to leading and first order 

We recover [vq-v>]^=0 and obtain d z V\-v + Vr • Vo = from equation (|3.18[) to 
orders —1 and to zeroth order respectively. From equation (|3.20[) we obtain 
3 z Mq = and ur = vo ■ v to order —2 and to order —1 respectively. To order —2 
equation (|3.19[) gives [t> ]+ = 0. 

To order —3, we have from (|4.22|) 

d t f%-u = 0, 

where 

Jl n * = -M T K{\ \d z <S> \ 2 + W(<f> ))d z Q v. 

This implies that J™\-v is constant in z. Furthermore, for any r such that 
r ■ v = 0, we have J™\ ■ r = 0. Hence = by gj25J| and this implies d z Q a = 0. 
Equation (|3.26[) gives to order —1 

Q = -d z (Ka(Q )d z $ ) + Ka(Q )W ($ )- 

Since d z Qo = we obtain = — d zz <&o + W'($o) again, which gives the profile for 
$o and the equipartition of energy (|4.9[) . Hence, we obtain the same solvability 
condition for $i from equation (|3.26[) : 

2^ = a(q ) K +[(G, l (c^(qo))-q c^(qo)))}l 

As previously, equation (|3.19[) then gives to order —1 

lp } 2 1 is-2l V ^D(v )} 2 1 i, = Ka{qo)v + Vrv(qo)- 

To order —2, we have from (14.221) 

d z J^-u = d z (J^-u)=0, 

where, thanks to <9 2 Qo = 0, 

j! n 1 t = -Mr(c r (g ))2W($ )(Vrgo + 5 2 giiv). 

This implies that 

d z (M r 2KW($ )d z Qi)=0. 
Integrating from — I to z and matching (|4.3p applied to $o gives that 

d z Qi = whenever |$ | < 1- 

4.5.3 Inner equations and solutions to second order 

To order —1, equation (|4.22p gives 

2 KW(* ) {d°c r (Q Q ) + Vo • V r c r (g )) + V 1 ■ ud z (2 KW(<5> Q )c r (Q )) 
where, using the already obtained results, = -M r {c r (Q ))2KW($ )V r Qo- 
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Proceeding as in Section 14.3.31 the left hand side yields 
d:(c r (q Q )) + c r (qo)Vr-v . 
For the right hand side, the integration from —I to +1 gives 

-VrYy^^-J^-H-i. 

where 

-V r • ^ 1^ J^j = V r ■ (MrVrgo) 

and m^gjl give 

-J^- l /|+; = -jr k ^|! = -(-Mi 2 )Vgo+Mi 1 )V g o)^=[jS]^. 
Hence we obtain the surface surfactant equation 

d* (c F (q )) + cf(q ) V r • «o = V r ■ (M r Vr9o) + [J%]\^- 

5 Numerical experiments 

In this section we report on numerical experiments that serve to support the 
above asymptotic analysis and illustrate that the proposed phase field models 
are able to describe phenomena that can be observed in physical experiments. 
Since the phase field approach to two-phase flow has been intensively studied 
already and the extension consists of accounting for the surfactant dynamics, 
the numerical experiments are designed to focus on the latter one. 

5.1 Surfactant adsorption dynamics in ID 

We first carefully investigate the adsorption of surfactants to interfaces in a 
one-dimensional setting where we exclude the effects of fluid transport (v = 0) 
and focus on the dynamics between bulk and interfacial surfactants. We assume 
that the surfactant is insoluble in fi^ and the sharp interface model is a variant 
of the Ward-Tordai problem defined on a bounded domain. For the phase field 
models we assume that tp is given, then the dimensionless equations of Model 
A simplifies down to (dropping the index * and the index 2 for the bulk phase) 

d t (e^c) - d x ( J-eo^c) = (y(c r ) - g'(c)) , 

d t (K5(v,d x p)c r ) - d x (J^S(<p,d x( p)d x c r ) - -/M(^,^)( 7 ' (c r ) - G'(c)) . 

For Model B we have one equation instead, 

dt (£(<p)c + K6(<p,d x (p)g(c)) -d x ^-^£,(tp)d x c+ ^-5(ip,d x tp)d x c) =0, 

and for Model C, we replace c,g(c),d x c by c(q),c T (q),d x q in the above equation. 
To support the asymptotic analysis we test 



35 



• the £-convergence of the profile of c(x,l); 

• the e-convergence of the profile of c r (0,i); 

• the £-convergence of |7'(c r ) — G'(c) | at x = 0,t=l. 

The third test only applies to Model A when j3 is scaled with e^ 1 , as the 
Dirichlct-type condition -f'(c T ) = G'(c) for instantaneous adsorption is enforced 
in the limit e — > 0. 

To measure the e-convergence of the profiles, we look at the difference 
\cpp — csi\ and |cp F — cF SI \ where c PF (x,t) and Cpf{x,£) are the interfacial and 
bulk densities of the phase field models respectively, while Cgj(i) and csi(x,t) de- 
note the interfacial and bulk densities of the sharp interface model respectively. 
We will be comparing {([5TT ]) , (53 j) } with Model A (a > 0) and {(EE]), (JO}} with 
Model A (a — > 0) and Model B. The numerical methods described in this section 
have been implemented using the software MATLAB, Version 7.11.0 (R2010b), 
[?]• 



5.1.1 Sharp interface model 



Set n = [0, 1] and T as the point x = 0, the dimensionless sharp interface model 
is 



d t c=^-d x 



in (0,1], 



d t c r = p^rd x c at x = 0, 



together with 



a a 
— o x c- 

Pe, 



(7'(c r )-G'(c)) at x = 



for non-instantaneous adsorption or 

c r (i)=.9(t) = (7')" 1 (G'(c)))atx = 



(5.1) 



(5.2) 



(5.3) 



for instantaneous adsorption. We impose the following initial-boundary condi- 
tions 

c(a: = l,f) = l, c(x,t = 0) = l, c r (t = 0) = c£. 

This is a version of the famous Ward-Tordai problem on a bounded in- 
terval, see [S3]. We solve the problem via a finite-difference scheme: Let 
= xi<---<2;Ar = l be a uniform discretisation of ft with mesh size h = l/N. 
Let At—l/Nf for integer NfGN be a time step and define t n = nAt for 
n = 0,...,Nf. Let 6 = At/(Pe c h 2 ) and denote c n (x) = c(x,t n ). Then given 
c™ = (c n (xi),...,c n (xN-i) 7 c n (xN)), the solution at time t n , we solved for c n+1 = 
{c n+1 (x 1 ),...,c n+1 {x N _ 1 ),c n+1 (x N )), which for {([EI]),^} satisfies 



/l + 26> -26 \ 

-0 1 + 26 -6 ... 





V o 



-6 1 + 26 -6 
1 / 



,n+i . 



/c"(x 1 ) + 2^( 7 '( C r '")-G'( C "(x 1 )))\ 
c n (x 2 ) 



c n (x N -i) 

1 
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and then 



„r,n+i _ J\n 



c 1 ■ n + 9h(c n+i (x 2 )-c n+i {x 1 )) 



For { (joTTjhfO)) }. we have to solve 

(eh -en o o \ 

-6 1 + 29 -0 ... 



... -6 1 + 29 -t 
V o 1 / 



,n+l 



/<?(c" +1 (^i))\ 



/<?(c"(*i))\ 
c"(x 2 ) 

c n (a;Ar_i) 

V i ) 



5.1.2 Phase field model 

We use the one-sided version for each of the above phase field models. We 
choose the potential W to be of double-obstacle type (hence K=-). This 
has the advantage that the phase field variable (p lies strictly in the interval 
[—1,1] and interfacial layer has constant width equal to £7r. The asymptotic 
analysis suggests that to leading order (p(x) =sin(|) for |ar| <£§ , and thanks to 
equipartition of energy S(f,d x ip) simplifies to 

%,M = {j > |0M( - )| 
The cutoff function £(*p(x)) is chosen to be 

^)=U(l + |(f)(3-(|) 2 ), 

U 

where y is the integer part of e%. 

For the discretisation we employ linear finite elements and the method of 
lines. Let At = ^- for integer Nf <E N be a time step and define t n = nAt for 
n = 0,...,Nf. Let Th be a uniform subdivision of the interval [—1,1] consisting 
of subintervals with size h. Let N be the number of vertices with coordinates 
denoted by {xi,--- ,xn}- Let Af be the set of vertex indices and for an index 
i£Af let LOi denote the neighbouring vertices connected to vertex i (i.e. Wi = 
{xi-i,Xi+i}). Furthermore, based on the functional form of 8 and £, we define 

Xh = {} € TV: there exists j € uii such that £(2^) > 0}, 
T^h ={i€Af: there exists j G such that <$(xj) > 0}. 

Let 

5 h :={w h eC ,0 ([-l,l]): u fc eP 1 ([x i) a; i+1 ]) ) i = l ) ... ) iV-l} 

be the discrete finite-element space. For t)GC°([—1,1]) wc define the interpola- 
tion operator IL h : C° ([-1,1]) S h to be 

N 

n fc fa):=5>(*i)Xi, 

i=l 



N<ef, 
\x\>e% . 



x>e* t 
M<ef, 

T < — F- 
— c 2 ' 
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where Xj( x ) denote the standard basis function such that \j £ 1, 1]) and 

Xj is a linear polynomial on each interval satisfying Xji x i) = $ji f° r an 

i,j = 1,...,N. Using the method of [T5], we can find the finite-element function 
c h ' n+1 (x) =c^(x,t n+ i) € S h such that c h ' n+1 {x.j) =0 if j £T>h and satisfying 

£ (/ 1 i n ' l ^" +1 ^)- 11 ' 1 ^"^)) + f ^i8)d xC r h n+1 d x X 3 
= -J 1 n h (l35( 1 i(c r h n )-Gi(c n h ))x J ), Vj€V h . 

The method for c£ +1 (x) = Ch{x,t n+ \) &S h is analogous, whereby c^ +1 (xj) = if 
j £ Xh and satisfies 

= £n h (J36{T'& n )-G'W))Xj), VjG^. 

For Model B, we seek c£ +1 G S h such that c^ +x {xj ) = if j £ X h U £> ft and satis- 
fying 

We remark that the scheme for Model C in this setting is structurally similar to 
the scheme of Model B. Hence in the subsequent one-dimensional experiments 
we will only implement the schemes for Models A and B, while Model C will 
be the subject of investigation in the two-dimensional experiments due to its 
two-sided nature. 

5.1.3 Numerics for Model A 

We observed the following regarding the choice of model parameters: 

• Interfacial Peclet number Per'. Fixing a = 0.2 and Pe c = l, we explored 
the effects of varying Per . For Per = 1 we observed that the profile for c r 
across the interfacial layer is linear when e = 0.2 or 0.1, but decreasing e to 
0.05 or 0.025 give a more uniform profile across the interface. Moreover, 
we can achieve a constant profile for larger values of e, i.e. e = 0.2 or 0.1, 
by decreasing Per to 0.01. 

• Bulk Peclet number Pe c : Fixing Per = 0.01, a = 0.2, we observe that 
the profile of c across the interface is linear for Pe c = 0.1 when e = 
0.2,0.1,0.05,0.025. When Pe c is increased to 10, we observe a constant 
profile in (— ef ,0) and a linear profile in (0,e-|). The size of these regions 
seems to be invariant for fixed Pe c as we reduced e from 0.2 to 0.025. 
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Model A Henry isotherm Model A Henry isotherm 




(c) (d) 



Figure 1: Model A ^-convergence for (a) the profile of c r (x = 0,t) and (b) the 
profile of c(x,t=l) with the Henry isotherm, (c) the profile of c T (x = 0,t) and 
(d) the profile of c{x,t — l) with the Langmuir isotherm. The parameter a is 
chosen to be 1. 
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Model A Henry isotherm, u -» Model A Henry isotherm, a -» 




(a) (b) 



Model A Langmuir isotherm, a -» Model A Langmuir isotherm, ra -» 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -1 -0.5 0.5 



(c) (d) 

Figure 2: Model A, e- convergence for (a) the profile of c r (x = 0,t) and (b) the 
profile of c(x,t=l) with the Henry isotherm, (c) the profile of c r (x = 0,t) and 
(d) the profile of c{x,t — l) with the Langmuir isotherm. The parameter a is 
chosen to be e. 
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These initial experiments with model parameters motivate the following 
choice for the convergence tests: We choose a=l, /?=-, Per = 0.01 and Pe c = 
10. The other parameters of the model are c\.j = l,c(a;,0) = l,c r (x,0)=0.05. 
The mesh size h is taken from {0.08,0.04,0.02,0.01,0.005} and the corresponding 
value of e is chosen from {0.4,0.2,0.1,0.05,0.025}. To ensure that the numerical 
scheme is stable, for each test we choose a time step At< h 2 . 

In the case of fixed a > we refer to Table [2] for the e-convergence in the 
difference in c r (0,l) and c(0,l) between the phase field model and the sharp 
interface model and Figure [1] for the e-convergence of the profiles. 



h 


e 




c PF (0,l)-(? SI (l) 




|c PF (0,l)-c SJ (0,l)| 


0.08 


0.4 


0.0974417 


0.0732749 


0.04 


0.2 


0.0419969 


0.0265120 


0.02 


0.1 


0.0163026 


0.0076752 


0.01 


0.05 


0.0058420 


0.0015298 


0.005 


0.025 


0.0022358 


0.0002207 


h 


e 




c PF (0,l)-c T SI (l) 




|cp F (0,1)-c sj (0,1)| 


0.08 


0.4 


0.0596860 


0.0963854 


0.04 


0.2 


0.0265857 


0.0364079 


0.02 


0.1 


0.0102234 


0.0115916 


0.01 


0.05 


0.0035830 


0.0030918 


0.005 


0.025 


0.0013697 


0.0009629 



Table 2: Convergence table for Model A, non- instantaneous adsorption (a = 1), 
Henry isotherm (top) and Langmuir isotherm (bottom). 

We also considered the scaling a~e (or p — e" 1 ) and from Figure [5] and 
Tables |3] and H] we observed the £-convergencc in the difference in c r (0,l) and 
c(0,l) between the phase field model and the sharp interface model. Further- 
more, we note that the maximum and mean difference of |7'(c r ) — G'(c) | in the 
interfacial layer decreases linearly as e — > 0. 

5.1.4 Numerics for Model B 

For Model B, since we have instantaneous adsorption, we can infer the difference 
of \cpp(0, 1) — cs/(0, 1)| from |cp F (0,l) — Cg 7 (l)| via the adsorption isotherms. 
Hence Table [5] displays only the difference |cp F (0,l) — Cgj(l)| for the Henry 
and Langmuir isotherms, in which we observe e-convergence along with Figure 
[51 The model parameters are chosen to be the same as in Model A. 

5.2 2D Simulations 

In this section we present some results of numerical simulations in two spatial 
dimensions in order to qualitatively illustrate the effectivity of our approach. In 
a first setting we expose a droplet of a fluid suspended in another fluid to a shear 
flow. Under moderate shear rates the droplet's shape attains a steady state. 
This shape changes in the presence of the surfactant. Of particular interest 
to us is the dependence of the shape on the isotherm. In a second setting we 
start with a droplet at rest (in particular, in equilibrium with respect to the 
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h 


E 








0.08 


0.4 


0.1191555 


0.1175129 


0.04 


0.2 


0.0685148 


0.0682569 


0.02 


0.1 


0.0383807 


0.0384228 


0.01 


0.05 


0.0209969 


0.0210621 


0.005 


0.025 


0.0114668 


0.0115106 


h 


e 


max| 7 '-G'| 


ave|V-G'| 


0.08 


0.4 


0.5882511 


0.1085532 


0.04 


0.2 


0.3540145 


0.0572062 


0.02 


0.1 


0.2061245 


0.0316161 


0.01 


0.05 


0.1128733 


0.0168467 


0.005 


0.025 


0.0594562 


0.0087458 



Table 3: Convergence table for Model A, instantaneous adsorption (a = e), 
Henry isotherm. 



h 


e 




c£ F (0,l) -<£,(!) 




|c PF (0,l)-c SJ (0,l)| 


0.08 


0.4 


0.0687143 


0.1452171 


0.04 


0.2 


0.0420765 


0.1452171 


0.02 


0.1 


0.0249919 


0.0506682 


0.01 


0.05 


0.0146093 


0.0292756 


0.005 


0.025 


0.0087232 


0.0173523 


h 


e 


max| 7 '-G'| 


ave|Y-G'| 


0.08 


0.4 


0.4014189 


0.0759004 


0.04 


0.2 


0.2347884 


0.0389953 


0.02 


0.1 


0.1326851 


0.0210856 


0.01 


0.05 


0.0711437 


0.0110897 


0.005 


0.025 


0.0370265 


0.0057192 



Table 4: Convergence table for Model A, instantaneous adsorption (a = e), 
Langmuir isotherm. 



h 


£ 


Henry 


Langmuir 


0.08 


0.4 


0.0938706 


0.0895642 


0.04 


0.2 


0.0616441 


0.0593439 


0.02 


0.1 


0.0336103 


0.0330060 


0.01 


0.05 


0.0172770 


0.0168309 


0.005 


0.025 


0.0083055 


0.0076996 



Tabic 5: Convergence table for Model B 




Figure 3: Model B e-convergence for (a) the profile of g(x = 0,t) and (b) the 
profile of c(x,t = l) with the Henry isotherm, (c) the profile of g(x = 0,t) and 
(d) the profile of c(x,t = l) with the Langmuir isotherm. 
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surfactant). Then we supply surfactant on one of the sides of the simulation 
box and investigate how far the droplet is sucked towards this side due to the 
Marangoni effect. As we are mainly interested in the effect of the surfactant on a 
qualitative basis we make convenient assumptions with respect to the two-phase 
flow, namely, that the fluids have the same mass densities and viscosities and 
that a Dirichlct boundary condition holds for the velocity. Also, the surfactant 
related parameters and data do not correspond to any specific species or systems. 

Both dynamic adsorption (Model A) and instantaneous adsorption (Model 
C) have been considered. In both cases, the Navier-Stokes-Cahn-Hilliard system 
was solved following the lines of [29] but we employed the double-obstacle po- 
tential for W(tp). The saddle point problem arising from (|3.29[) and (|3.33|) has 
been solved with a preconditioned GMRES [15] ■ For the phase field equation 
(|3.31[) together with (|3.32p in form of a variational inequality we have employed 
a Gauss-Seidel type iteration as described in [S]. 

For Model A, we always considered Fickian diffusion by setting M^l(c*) = 



l/(G'*'(c*)Pe c ,i) and M r ,»(c£) = l/(7^(c^)Pe r ,j). We also replaced 6*((p,V*<p) 



by 2W(<p)/e* in the surfactant equation (|3.35[) which effects the validity of the 

energy inequality but doesn't change the result of the asymptotic analysis. The 
reason is that the method developed in [TO] can directly be applied. We leave 
a careful study of the impact of the gradient term for future investigations. In 
analogy to [TO] a method for the degenerate bulk surfactant equations (|3.34[) 
has been developed. The methods have been implemented using the software 
ALBERTA, Version 2.0.1, gS]. 

In the surfactant equation (|3.37[) for Model C we assumed constant mo- 
bilities, Mc)* (c* (q*)) = 1/Pe c ,j and Mr,*(c» (q*)) = 1/Per, and we also replaced 
<5*(vj, V,^) by 2W(tp)/e* for not having to deal with V*c^ in the diffusion term. 
Whenever no closed formula for c£, c*\ or as a function of q* was available 
we employed a Newton method. In the same way we also dealt with the nonlin- 
ear system of equations emerging from the finite element discretisation of the 
surfactant equation. 

With regards to parameters and functions appearing in non-dimensional 
equations of the phase field models we have in both settings: K = 2/tt, X p = 1, 
Xr, = l, Ca=0.1, 



and ^(v) = 1 — £i(¥>) where we set ^((p) = if \<p\ > 1. 
5.2.1 Droplet in shear flow 

On the domain 0= [— 5,5] x [— 2,2] Cl 2 the velocity was initialised with 
v(xi,X2,0) = 0. On the upper and lower boundary {a;2 = 2} and {x2 — — 2} we 
then increased the velocity linearly in time to v(x\,X2,t) = (x2/2,0), t > 0.1. On 
the two sides {x\ — — 5} and {xi=5} we imposed the condition v{x-i,X2,t) = 0. 
The phase field was initialised with ip(x,Q) = i[i((\\x\\2 — where 





1<<P, 
-l<tp<l 

¥><-!, 




for z>\ , 
for |z|<f, 
for z<-§, 



(5.4) 
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Figure 4: Droplet in shear flow: Zero level sets of ip for several isotherms, 
e = 0.0565685425 f«0.08/v^, t=10. The right graph displays a zoom into the 
square indicated on the left graph. 




Figure 5: Droplet in shear flow: Interface surfactant density c£ (left) and surface 
tension cr»(c£) (right) plotted over the angle formed by the line from the centre 
to a boundary point and the x-axis for several isotherms, e — 0.0565685425 « 
0.08/V2, i = 10. 
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angle angle 



Figure 6: Droplet in shear flow: Surface tension er*(c£) at the tips of the droplet 
(left) and difference of surface and bulk chemical potentials 7*(c*) — G^(c^) 
(right) plotted over the angle formed by the line from the centre to a boundary 
point and the x-axis for several values of a*, £ = 0.0565685425 0.08/ y/2, £=10. 

which yields a circular diffuse interface of radius one and centre m = (0,0). 
Furthermore, we set Re = 0.1 and m*((p) = |(1 — f 2 )+- 

We investigated Model A with Per = 2.5, Pe c j = 2.5, and a* = 1 for i = 1,2 
for the following isotherms, see Table [T] (assuming the same free energies in the 
two bulk phases, thus dropping the index): 

• Langmuir (5 = 0.2, <r = l, if =10); 

• Frumkin (5 = 0.2, a = l,K=W, A = 0.4); 

• Frcundlich (B = 0.2, a = l, A' = 10, A = 1.5, A c = 1.0). 

The initial bulk surfactant density was c* =c^ = l/(10e)~ 0.03679, and the 
interfacial surfactant density was the equilibrium value (thus, depending on 
the isotherm). 

At time t = 10 the droplets seemed to have attained stationary shapes. These 
are displayed in Figure 0] for several isotherms. For our parameters we found 
that the Langmuir isotherm leads to the least deformed shape while the shape 
for the Freundlich isotherm is most deformed when comparing with the initial 
circular shape. A common measure for the deformation is the Taylor deforma- 
tion parameter DTay = {L — B)/(L + B) where L and B are the maximum and 
the minimum distance to the centre, respectively. We obtained the following 
values: 



isotherm 


Langmuir 


Frumkin 


Freundlich 


D T ay 


0.143298 


0.148370 


0.160821 



In Figure [5] we display the surface surfactant density and the surface tension 
along the interface between the two fluids which qualitatively reveal the usual 
distribution, for instance, compare with [36]. 

We also investigated a change in the adsorption parameter (both always 
equal for the two phases, whence we drop the upper index). The impact on the 
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Figure 7: Marangoni effect on a surfactant laden droplet due to the provi- 
sion of surfactant at the boundary. Computed fields 2W ((p)c^ (q„) (left) and 
£1 (<fi)c^ ((/*) (right) are plotted over the domain O = [—3,3] x [—2,2] (x-axis from 
left to right, y-axis from front to rear, z-axis or height indicates the value of 
the field) at times t = 0,10, 40, 100 (top down) for a simulation performed with 
the Frumkin isotherm data (see Section [5.2. 2[) and e = 0.12. The data range is 
between 0.0 (blue) and about 0.585 (red). 

shape is small in comparison with the isotherm. For the Langmuir isotherm, we 
obtained the deformation parameters 



adsorption parameter 


a* = 2.0 


a* = 1.0 


a* = 0.5 


Dray 


0.143395 


0.143298 


0.143241 



In Figure[B]thc difference of the chemical potentials at the interface is displayed, 
revealing the expected convergence to zero when the adsorption parameter a* 
decreases. 

5.2.2 Marangoni effect 

We now consider the domain fi = [— 3,3] x [— 2,2]. Both velocity and pressure 
are initialised with 0, and this is also the Dirichlct boundary condition for the 
velocity. For the phase field we set <^(x,0) =ip((\\x — m\\2 — l)/^) with if; given 
as in (|5.4|) and m= (0.5,0) which corresponds to a circular diffuse interface of 
radius one around m. The Reynolds number is Rc= 10 and we chose m*(<p) = 

(1-^)+. 

Simulations were performed with Model C where we set Pe C; ; = Per = 10.0, 
i = 1, 2 and used the following free energies (again, the free energies in the two 
bulk phases are assumed to be the same so that the index is dropped): 
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• Langmuir (B = l, o-q = 2, AT = 2.5); 

• Frumkin (B = l, a = 2, K = 2.5, A = 0.4); 

• Frcundlich (B = l, a a = 2, K = l, N = 1.5, A c = 0.6). 

The field q* was initialised such that c* (<?*) ~c^\q r ) = 0.1 . During the time 
interval [0,0.1] we linearly increased g* on the boundary {x\ = — 3} such that, 
at t = 0A, ci 1) (^) = 0.5. 

As a consequence, the droplet moved in —x\ direction towards the source 
of the surfactant as exemplary illustrated in Figure [7] for the Frumkin isotherm 
data. Initially at rest, the supply of surfactant on the boundary leads to a 
surfactant gradient at the interface of the droplet. Since <r* is decreasing in 
the related Marangoni force Vrc*(c£) points into the opposite direction and, 
thus, leads to a drift towards the source of the surfactant. In the long term, the 
system reaches a steady state again with spatially homogeneous distributions 
of the surfactant in both phases and on the interface, which is fairly achieved 
at time t= 100.0. For our choice of parameters the Frcundlich isotherm lead 
to the most significant displacement d xi along the x\ axis while the Langmuir 
isotherm lead to the least significant displacement: 





Langmuir 


Frumkin 


Freundlich 


d Xl 


-1.055512 


-1.087783 


-1.114869 



6 Appendix 

We use the following result from Alt [3] to reformulate the strong form of the 
surfactant equations (|2.8|) . (|2.11[) . (|2.12p into an equivalent distributional form. 
Let D'(f2) denote the space of distributions on f2. 

Theorem 6.1 (Alt [3J Section 2.7 &: Theorem 2.8) Given an open set 
D C R x M. d consisting of two open sets H*- 1 -* and f^ 2 ) separated by a smooth 
evolving hypersurface T, in particular, TdD has no boundary within D. For 
(t,x)eT we let Vi{t,x)e {T x (T(t)))^ dW l be the external unit normal ofQ (l) (t). 
Then V\ + 1/2 = 0. Denote by Xaw >Xfi< 2 ) i^r the following distributions: 

[ fdxnm= [ [ /(*,*)> [ fdfir= [ [ f(t,x). 

Then a single balance law is an equality of the form 

d t E + V-Q = F inV'(D) (6.1) 

with distributions given by 

i=1.2 i=l,2 £=1,2 

where &^\qf ,f^ ■ fi® ~ >K and e r ,qj,/ r :T— >R are smooth functions. Then 
the distributional law (|6.1[) is equivalent to the following: 

1. For i = 1,2 in ftW; 
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2. For all (t,x)eT: 

(q r -e r u r )(t,x)eT x (T(t)). 

3. OnT: 

d t e v + u T -Ve T -e r K V -u r +V T -{q T -e r u T ) = f + ^ (gW -e®u r )-Vi, 

i=l,2 

where Ur is the unique velocity vector such that 

T (t)X) r=span{(l,u r {t,x))}®({0}xT x T(t)), 
and Kp is the curvature vector defined by 

Vr ■ n = —Kr-n, 

for spatial normal vector fields n(t,x)G (T x T(t))' L . 

For the reformulation, we assume as in [48] that c r is extended off T constant 
in the normal direction, hence Vrc r = Vc r . Define 

ji = ^(7'(c r )-Gi(c (1) )), , 2 = -L( 7 '( c r )-G 2 (c( 2 >)), 

then by the definition of d'(-), the divergence- free property of v and that 
V7'(c r )=7"(c r )Vc = 7"(c r )V r c r = Vr7'(c r ), equation (gUP can be written 
as 

9 t c r + V r -(c r - i ;-MrV7'(c r )) = -(ii+j 2 ). 

Choosing e W = gj i) = /W = for i = l,2, l<j<d and e r = c r ,q r = c r v~- 
Mr'Vj' (c r ), f r = — {jx+ji)- Theorem 16.11 implies that the distributional form 

d t (S r c r ) + V • (S r c r v - A/ r 5rV7'( C r )) = -5 r (ji +j 2 ) (6.2) 

is equivalent to 

d t c r +u r ■ Vc r - c r n r ■ u r + V r ■ (c r v - il/ r V r 7'(c r ) - c r u r ) =-(ji +32) on V. 

We have Vr-(c r «r) = -c r Kr'Mr and ur = {v-i>i)vi implies v = up + v T . Fur- 
thermore, Vr • (c T v) = Vrc r ■ v T + c r Vr ■ i>. Hence equation (|2.16[) is equivalent 
to (J2JTJ). For i = 1, choose e (2) = gj 2) = = / (2) = e r = = for 1 < j < d and 
e (1) =c (1) ,/ r =ji,g (1) =c (1) i)-Ad 1) VGi(c( 1 )). Then the distributional form 

dt (xnw c (1) ) + V • (xnci) c (1) u - Xn w VGi (c«)) = frji (6.3) 

is equivalent to 

5*(c (1) ) + V-(cWt;-MPVG' 1 (c (1) )) = 0, in ftW, 

•z/i =j'i, on r. 



49 



Similarly, choosing = = /M = = e r = = for 1 < j < d and = 
c (2) 5 f r =j 2 , q( 2 ) = c^v - M c (2) VG" 2 (c( 2 )). Then the distributional form 

ft (Xnw c (2) ) + V • (x (2, « - xnw M i 2) VGf, (c< a > )) = <5 r . ?2 (6.4) 
is equivalent to 

d t {c^) + W-{c^v-M^G' 2 {c^))^Q, in 

-MWvG , 2 (cW)'V 1 =j 2 , ouT 

as V2 = —v\ . Thus the bulk and interfacial surfactant equations can be refor- 
mulated in the distributional forms f|6 . 2[) — (|6.4[) . 
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